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1.  INTRODUCTION 


Aeronautical  Research  Associates  of  Princeton,  Inc.  (A.R.A.P.) 
has  been  developing  a program  for  computing  the  detailed  low- 
level  atmospheric  distributions  of  velocity,  temperature,  moisture, 
refractive  index,  and  the  turbulent  variances  of  these  quantities 
for  marine  environments.  Our  approach  to  this  problem  has  been 
to  use  the  invariant  second-order  closure  model  of  turbulence 
developed  by  Dr.  Coleman  duP.  Donaldson  and  his  associates  at 
A.R.A.P.  over  the  past  few  years.  The  fundamentals  of  this  ap- 
proach are  given  by  Donaldson  (1973)  . A review  of  the  status  of 
this  model  as  applied  to  a wide  variety  of  turbulent  flow  problems 
is  given  by  Lewellen  (1977) . Particular  applications  of  the  model 
as  applied  to  atmospheric  problems,  including  comparisons  with 
experimental  data,  are  documented  in  Refs.  2-6. 

Reference  7 is  a technical  report  detailing  the  model  devel- 
opment, sample  calculations,  and  verification  comparisons  made 
under  our  initial  NASC  contract.  It  describes  the  addition  of 
humidity  and  the  second-order  turbulent  correlations  involving 
htimidity  as  variables  to  our  dry  atmospheric  boundary  layer  program 
previously  .developed  for  the  Environmental  Protection  Agency. 

Using  the  predicted  distributions  of  temperature,  humidity,  and 
pressure,  a calculation  of  the  modified  refractive  index,  M , 
has  been  incorporated  in  the  program.  Local  minimums  in  the  M 
distribution  with  respect  to  altitude  directly  indicate  the  presence 
of  a radar  duct.  Since  we  are  predicting  the  second-order  correla- 
tions between  the  turbulent  fluctuations  in  temperature  and 
humidity  as  well  as  the  average  scale  of  the  turbulence,  we  have 
available  the  information  to  also  compute  the  structure  of  the  fluc- 
tuations in  refractivity . 

The  second  year  of  the  research  called  for  two  major  modifica- 
tions to  the  model  described  in  Ref.  7.  These  were  to  (a)  increase 
the  dimensions  of  the  program  to  a two-dimensional,  unsteady  cal- 
culation to  permit  the  prediction  of  shoreline  conditions  de- 
veloping in  time,  and  (b)  to  incorporate  the  radiative  flux  diver- 
gence term  into  the  one-dimensional  system  of  equations  in  a 
coupled  manner.  These  developments  are  detailed  in  Px.ef.  3. 

The  two-dimensional,  unsteady  version  of  the  model  has  been 
used  to  calculate  the  typical  variation  in  the  coastal  planetary 
boundary  layer  (Ref.  9).  The  resulting  diurnal  variation  in  the 
sea  breeze  induced  by  the  strong  stability  differences  in  the 
boundary  layer  response  over  the  land  and  that  over  the  water 
produces  a strong  as3n3imetry  between  the  sea-breeze  and  the  land- 
breeze  circulation  patterns.  In  previous  sea-breeze  models,  it 
was  necessary  to  impose  eddy  viscosities  which  were  a strong 
function  of  time  and  space  to  gain  this  as3rmmetry.  In  the  present 
model,  it  is  achieved  without  the  need  to  introduce  any  new  em- 
pirical information. 

The  last  year's  effort  has  been  divided  between  exercising 
lodel  for  verification  purposes,  and  extending  its  capability. 
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Extensions  to  the  model's  capability  include:  (1)  the  removal  of 
the  hydrostatic  approximation;  (2)  the  incorporation  of  condensa- 
tion and  radiation  into  the  two-dimensional,  unsteady  model; 

(3)  the  incorporation  of  terrain  variation  into  the  surface 
boundary  conditions ; and  (4)  the  ability  to  determine  the  influence 
of  particle  size  on  the  turbulent  transport  of  aerosols.  These 
developments  are  discussed  in  Chapter  2.  Some  of  the  sample 
calculations  are  detailed  in  Appendix  A,  a paper  of  our  modeling 
results  which  has  been  accepted  for  publication  in  the  Journal  of 
the  Atmospheric  Sciences.  These  sample  calculations  demonstrated 
the  strong  dynamic  interaction  between  turbulent  transport  and 
thermal  radiation.  Other  sample  calculations  are  discussed  in 
Chapter  3 and  Appendix  B. 

We  have  also  investigated,  analytically,  some  of  the  relation- 
ships between  parameters  governing  distinct  fronts  in  the  planetary 
boundary  layer,  as  an  aid  toward  the  detailed  computation  of  the 
structure  of  such  fronts.  This  analysis  is  discussed  in  Appendix  C. 


i 
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2 . MODEL  EXTENSIONS 


Several  important  modifications  to  the  model  have  been  made 
during  the  last  year.  We  began  with  the  removal  of  the  hydrostatic 
approximation  from  the  two-dimensional  program.  This  was  ac- 
complished by  eliminating  the  pressure  as  a primary  variable  and 
formulating  the  equations  in  terms  of  streamfunction  and  vorticity. 
The  numerical  implementation  of  this  modification  is  not  described 
in  the  present  report,  but  the  results  of  a computation  where  the 
influence  of  nonhydrostatics  is  important  is  described  in  Section 
3.2.  The  same  two-dimensional  calculation  includes  the  effect  of 
condensation  and  evaporation.  This  involved  a straightforward 
carry-over  of  the  same  change  of  phase  model  as  described  in 
Appendix  A. and  applied  to  the  one- dimensional  code.  Another  mod- 
ification to  the  numerics  not  described  in  detail  herein  was  to 
provide  a stronger  coupling  between  the  turbulent  correlations. 

This  permitted  the  simulated  time  step  in  our  one-dimensional  code 
to  be  increased  by  as  much  as  a factor  of  5 . Model  extensions 
which  have  been  formulated  but  are  not  represented  in  any  of  the 
sample  calculations  are  discussed  in  the  following  sections. 

2 . 1 Implementation  of  Atmospheric  Radiation  Transport  in  the  Two- 
Dimensional  Model 


The  A.R.A.P.  atmospheric  radiation  mode],  described  in  Lewellen, 
et  al. , 1976  has  been  implemented  with  the  two-dimensional  tur- 
bulence  models  for  the  situations  where  vertical  radiative  trans- 
port is  dominant.  This  includes  most  atmospheric  boundary  layer 
motions  including  those  with  clouds  consisting  of  stratus  decks, 
closely  packed  cumulus  and  fog  layers.  Typically  we  require 
(1)  that  cloud-cloud  horizontal  radiative  transport  be  negligible 
compared  to  vertical  radiative  transport  or  horizontal  advective 
and  turbulent  transport,  and  (2)  that  radiation  from  the  sides 
of  cumulus  cells  to  the  ground  or  to  space  be  negligible  compared 
to  vertical  radiative  transport  in  the  cumulus.  If  we  denote  the 
average  depth  of  a cloud  as  5 and  its  characteristic  breadth 
as  L , then  a condition  for  vertical  dominance  of  radiative 
transport  may  be  formulated  by  requiring  that  radiative  cooling 
to  space  through  the  cloud  top  be  much  greater  than  that  through 
the  sides.  In  terms  of  the  fractional  cloudiness  R , this  con- 
dition is 


(6/L)  ^ - 1) 

(6/L)2  1)2 
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(2.1) 


When  the  clouds  are  widely  spaced  (R  <<  1)  this  condition  just 
requires  that  the  "isolated"  cloud  have  a small  aspect  ratio 


6/L  <<  1 


(2.2) 


On  the  other  hand,  as  the  sky  is  more  filled  with  clouds  (R  -»■  1) 
the  aspect  ratio  condition  is  relaxed  due  to  the  blocking  of 
radiation  to  space  by  adjacent  clouds. 


Within  this  approximation  the  vertical  radiation  field  at 
a local  horizontal  position  is  determined  by  the  distribution  of 
fluid  properties  in  the  vertical  at  that  horizontal  position. 

The  horizontal  variation  of  the  radiation  field  thus  is  a function 
of  the  horizontal  variation  of  fluid  absorption,  scattering  and 
emission  properties. 

It  is  interesting  to  note  that  the  condition  (2.1)  is  essen- 
tially a requirement  that  radiative  loss  (or  gain)  in  the  horizontal 
direction  be  small  compared  to  the  vertical.  As  such  it  is  also 
the  appropriate  condition  for  validity  of  the  vertically  struc- 
tured direct  solar  radiation  scattering  and  absorption  model  in 
which  only  vertical  back- scattering  of  incident  solar  energy  by 
water  droplets  and  aerosols  is  included.  The  required  condition 
for  validity  of  the  vertical  radiation  model  in  a two-dimensional 
calculation  is  then  that  scattering  of  solar  radiation  off  the 
sides  of  clouds  be  negligible  compared  to  cloud  top  scattering. 

2 . 2 Two-Dimensional  Computations  with  Terrain 

We  now  describe  the  procedures  for  the  treatment  of  non- 
uniform  terrain  in  the  A.R.A.P.  two-dimensional  atmospheric 
boundary  layer  calculations.  We  consider  a rectangular  y-z 
domain  as  shown  in  Fig.  2.1  in  which  a rectangular  grid  of  m x n 
points  exists.  A nonuniform  terrain  boundary  3T  cuts  through 
this  grid  which  we  shall,  for  simplicity,  assume  passes  through 
existing  grid  points  as  shown  in  Fig.  2.1.  (The  extension  to 
terrain  which  does  not  fall  exactly  on  the  existing  grid  points 
is  straightf orchard  and  will  not  be  discussed  here.)  We  let 
z ■ Zm  be  the  coordinate  of  the  terrain  and  we  let  Zq  be  the 
roughness  height  measured  from  Zjjj  . Let  N be  the  total  num- 
ber of  interior  points  of  the  domain  which  lie  on  the  terrain 
boundary  3T  . 

The  fundamental  mean  variables  in  two-dimensional  calculations 
consist  of  the  vorticity  r\  , potential  temperature  6 , and  the 

streamf unction  4/  describing  the  mean  velocity  field.  These 
variables  are  supplemented  with  the  turbulent  correlations  and 
their  corresponding  equations.  The  streamfunction  and  vorticity 
are  defined  on  the  full  rectangular  grid  but  all  other  variables 
are  defined  on  only  the  grid  points  where  fluid  exists  outside 
the  terrain  and  on  the  boundary.  Our  purpose  in  defining  the 
streamfunction  ip  and  vorticity  n on  the  full  grid  stems  from 
the  need  to  rapidly  solve  the  streamfunction  equation 

- -n  (2.3) 

at  each  time  step.  For  this  purpose,  we  utilize  a rapid  direct 
cyclic  reduction  algorithm  which  requires  a uniform  boundary 
domain.  Thus,  we  solve  Eq . (2.3)  on  the  full  rectangle  including 
the  region  lying  "inside"  the  terrain.  This  procedure  will  not 
invalidate  our  solution  for  the  flow  outside  the  terrain  because 
data  inside  the  terrain  will  not  be  utilized  and  the  solution 
will  be  constructed  such  that  the  correct  terrain  surface 


conditions  on  tp  in  particular  and  all  other  variables  in  gen- 
eral will  be  impressed  on  the  embedded  terrain  points  by  the 
method  we  now  describe  (Roache,  1972). 


The  boundary  conditions  on  the  rectangle  are  specified  in 
the  usual  fashion  for  all  variables  except  on  that  portion  of 
the  rectangle  lying  within  the  terrain.  On  that  portion,  no 
boundary  values  are  required  except  for  the  streamf unction  ip  . 

The  value  to  be  imposed  on  p is  arbitrary  on  that  bovindary 
segment  and  we  simply  make  it  a no-flow  boundary  arbitrarily  setting 
ij;  = constant  =>  . Inside  the  terrain  we  set  arbitrary  values 

of  the  vorticity  n . typically  n = 0 . 


Let  us  first  consider  the  method  for  determining  the  stream- 
function  which  satisfies  Eq.  (2.3)  within  the  domain  of  the 
fluid  outside  the  terrain  and  the  boundary  condition  of  no-flow 
through  the  terrain,  i.e.. 


constant  = on 


3T 


We  represent  the  streamfimction  (i  = 1,2, 


m X n ) as 


p . ^ p + y a . G . . 


N 


where 


satisfies 


(0) 


(2.4) 


(2.5) 


aj  is  an  N length  column  vector,  and  Gij  (i  = 1,2,...  m x n ; 
j = 1,2,...  N)  is  an  (m  x n)  x N matrix  which  satisfies 


(2.6) 


where  5ij  is  the  Kronecker  delta  and  the  index  j ranges  over 
all  the  points  j = 1,2,...  N which  lie  on  the  terrain  boundary. 
The  matrix  function  Gij  is  a discrete  Green's  function  for  the 
boundary  and  is  a discrete  difference  operator  approximating 

the  Laplacian  operator. 


Equations  (2.5)  and  (2.6)  are  solved  on  the  full  rectangular 
grid  with  boundary  conditions  on  ij;(0)  given  by  those  required 
by  the  physics  of  the  problem  outside  the  -terrain  and  by 

on  the  portion  of  the  rectangle  lying  inside  the  terrain. 
The  Green's  functions  Gij  satisfy  either  Gij  =0  or 
3G^j/3n  * 0 on  the  rectangular  boundary,  depending  upon  whether 
p takes  a Dirichlet  or  Neumann  condition.  N + 1 executions 
are  required  for  the  cyclic  reduction  algorithm;  One  for  p^^^ 
and  N for  the  G^j  . 
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The  constants  aj  are  determined  by  requiring  that 
take  on  the  specif ied"^ value  on  3T  . Thus  on  the  N points 

lying  on  9T  we  require 

/ns  N 

i = 1,2,  . . ..'J 

where  G^j  is  a portion  of  matrix  Gij  only  containing  contribu- 
tions from  Gij  for  i corresponding"^  to  points  lying  on  9T  . 

The  matrix  G^j  is  thus  an  N x N matrix.  The  constants  aj 
are  deternined"*as 


+ .1 


a.  = - I g:1  - O (2.8) 

J i=l  ^ ^ 

We  note  that  the  matrices  Giji  , Gij  , and  Gl}  may  be  computed 
once  for  all  time  since  they  depend  only  upon  the  geometry  of 
the  domain  (including,  of  course,  the  terrain  boundary). 

After  each  time  advance  in  which  an  updated  n field  is 
obtained,  we  solve  Eq.  (2.5)  for  . We  then  determine  the 

solution  satisfying  the  terrain  boundary  condition  and  the 

remaining  boundary  conditions  on  the  rectangular  region  as 


'^i 


N N 


'^i 


I I 0-  (^^ 

j = l i=l  ^ 


(0) 


- V^ij 
«-l 


utilizing  the  known  and  stored  values  of  G.^  and 


(2.9) 


The  vorticity,  temperature  and  turbulence  variables  are 
solved  only  in  the  region  outside  the  terrain  utilizing  an  al- 
terna ting -direction- imp lie it  technique  which  does  not  require 
uniform  boundaries.  Boundary  conditions  on  the  rectangular 
portion  are  conventional,  but  boundary  conditions  on  the  terrain 
boundary  9T  require  special  treatment.  As  on  other  no-slip 
boundaries,  we  do  not  integrate  the  equations  all  the  way  to  the 
wall  (which  would  require  exceedingly  fine  and  wasteful  mesh 
spacing).  Instead,  we  require  that  the  variables  take  on  their 
constant  flux  region  values.  For  the  vorticity  this  requirement 
results  in  a linear  relationship  between  the  boundary  value  of 
n and  ip  on  the  first  interior  point.  For  the  turbulent  cor- 
relations, the  requirements  is  that  their  normal  gradients 
vanish. 


Preliminary  results  demonstrate  the  feasibility  of  the  Green's 
function  approach  to  give  the  correct  streamfunctions  throughout 
the  computational  domain.  The  procedure  yields  a fairly  realistic 
contour  pattern  for  the  streamfunction  at  flow  separation  over  a 
two-dimensional  mountain,  as  shown  in  Fig.  2.2.  Here  the  flow 
is  entering  from  the  left,  separating  down  the  slope,  and  reattaching 
further  downstream.  An  added  divergence  of  the  flow  at  the  top 
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of  the  calculation  domain  enhances  separation.  The  separation 
point  along  the  slope  is  fairly  v;ell  held,  but  the  bubble  itself 
seems  to  push  the  reattachment  point  slowly  downstream.  Some 
questions  must  still  be  answered  with  regard  to  the  best  way  to 
apply  the  vorticity  and  turbulence  boundary  conditions  along  the 
interior  boundary. 


2 . 3 Turbulent  Transport  of  Aerosols 

In  previous  reports  (Donaldson,  1973,  Lewellen  and  Teske, 
1975)  second-order  closure  has  been  applied  to  gaseous  diffusion. 
The  purpose  here  is  to  apply  the  same  technique  to  the  diffusion 
of  finite  size  aerosol  particles,  in  order  to  determine  what 
conditions  influence  their  turbulent  diffusion  and  how  the  dif- 
fusion of  such  particles  may  be  computed. 


Following  Marble  (1970)  the  equations  of  motion  for  an  in- 
compressible flow  containing  a cloud  of  particles  of  concentration  c 
may  be  written  as 


3u. 

377=0  (2.10) 


3c 

Jt 


+ 


S 


(2.11) 


9u 

Tt 


3u.u.  T ^ 

1 = i ^ 


p 3x.  S K 


J 


+ F. 
1 


V 

3cua^  3(cuA^UAj)  _ 


3t 


3x . 


= -c 


+ cF . 


(2.13) 


The  characteristic  time  is  the  time  required  by  a particle 

to  reduce  its  velocity  relative  to  the  gas  by  e"-^  of  its 
original  state,  in  the  unaccelerated  state.  For  a significant 
range  of  particle  radii  may  be  described  by  Stokes  flow 


with  the  density  of  the  aerosol  particle,  and  a its  radius. 

We  will  follow  the  same  procedure  of  Donaldson  (1973)  to 
obtain  equations  for  the  turbulent  transport,  u^  c ' , of  the 

aerosol.  When  the  variables  are  decomposed  into_in  ensemble 

mean  and  a fluctuating  contribution,  i.e.,  c=c+c'  , 

u.  * u.  + u'  , etc.,  then  Eqs.  (2.10)  to  (2.13)  may  be  averaged 

Ai  A^  A^ 
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k + J 1 

at  ax. 


J -J  = _ 12  * 11  ^ ^ ^1  I-  — CU.  * ' 

axj  “a, “A.  3Xj  -Aj  3Xj 

*“a,  3c'u'  ^ (PTiJ-JMT) 

A„  ~ iw  ~ X 


- C Ua  ^ 

A.  ax. 

+ 5 


(2.20) 


An  equation  for  the  correlation  ul  ul  appearing  in  Eq, 

Ai  Aj 

(2.20)  may  be  obtained  by  manipulating  the  equation  for  u^ 

"i 

derived  by  subtracting  Eq.  (2.19)  from  Eq.  (2.18). 


Ai  A. 


^’^A.^A.  “"^A.  ““A, 

^A^  ax^  ^ “ ■ ^Aj’^Aj^  axJJ 


A ^ ■ \ ^ 


- "A,  ■ 


j 3- 


1 J 


- “Aj 


(2.21) 


An  equation  for  the  cross  correlation  u'  u!  may  be  obtained 
from  Eqs.  (2.12),  (2.16),  (2.18)  and  (2.19)  to  be 


aui  u! 

Ai  J 


at  \ ax^  ■ ^A^'^  ax^ 

3Ua.  ^“a. 

■ ^ Nr  ■ rf 

- \ lij  + ^ ^ 


+ V terms 


(2.22) 
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An  equation  for  c'u^  can  also  be  generated  from  Eqs . (2,11), 
(2.15),  (2.12)  and  (2.19)  to  be 

3c ‘ul  _ 3c 'u!  _ _ r ~- 

—5-= — + u.  —5 — - = -(u.  - u.)  ul  ^ uTuT 

3t  j 3x.  'A.  J 1.  3x.  1 A.  3x. 

•J  J ^ w 

3u.  3c ' ulul  ^ 

" ^'^'2  ^ 3x^  ^ * ^^A.  ■ 3x7 

«J  w J 

- ='  ^ (Saj^  - ij.) 


+ — (c ' ‘‘u'  - c'  ) + V terms 


(2.23) 


The  extent  of  the  new  modeling  required  to  close  the  system 
by  carrying  the  full  set  of  equations  for  c 'uI  , c'u:  , 

uT~uT  , ul  ul  , along  with  uTuT  may  be  assessed  from  Eqs. 

"i  "j  J 1 j 

(2.20)  through  (2.23).  The  list  is  quite  extensive: 


'^A.^A.  3x,  ’ ^'^A.''^j^'\3x.  ’ 3x,  ’ '^aX^  ’ 

r j K j j K 1 


(“1,  - “i)  “j  u:||- 

IX  .L  .L  .L  J 

along  with  the  terms  involving  S'  and  v . Rather  than  undertake 
this  modeling  task,  we  will  take  advantage  of  the  fact  that 
Tv  associated  with  most  aerosols  of  interest  is  much  shorter 
than  the  characteristic  time  of  the  atmospheric  turbulence, 

A/q  . We  will  assume  that  u^  may  be  expanded  as 


Ua_  “ + — (^^1.  ) + 0(T^q/A)' 


(2.24) 


If  this  expansion  is  carried  through  Eqs.  (2.20)  through  (2,23), 
the  leading  order  terms  reduce  to 
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Equations  (2.25)  through  (2.23)  must  also  be  consistent  with 
the  equation  for  uj^uj  computed  from  Eqs.  (2.12)  and  (2.16). 


3u|ul  _ 3uTuT 

— — ^ + u,  1 = -uTu 

3t  K ax,  1 1 


. \1  —I  - u . Ui  -■ 


axj, 


- “1  Hj  * - “j)  + 

+ U^^uj)  + - i^)  + 


+ V terms 


(2.29) 


First,  let  us  consider  the  case  where  c <<  1 and  the  fluid 
turbulence  may  be  decoupled  from  the  particle  motion.  In  this 
limit  the  terms  involving  c and  c’  may  be  dropped  from  Eq. 
(2.29).  Let  us  also  assiome  S = 0 , and  note  that  u^  - u. 

“ • Then  subtracting  Eq.  (2.29)  from  (2.26)  ind  (2.27) 

leads  to 

“ - “i  Hj  - p <2;^^  - - np:j^)  + c (2.30) 

° ^ ‘ “a.  X ^ 

Also,  from  Eq.  (2.24)  it  is  possible  to  show  that 

^ + 0(T  q/A)^  (2.32) 

i-J  i j ij  J 

Thus,  Eq.  (2.30)  leads  directly  to 

- 3T3T  + T^Cui  111  + O 

J 3 

I^Jhichwith  our  usual  model  for  u|3p ' /3Xj  and  e gives 

“a^Aj  - ^ ^ <2.33) 

Similarly,  from  Eqs . (2.25)  and  (2.28) 


C * U J = c ' i_i  1 + T c'  

IV  3x^ 


c Uj;  (1  + AT.^q/A) 


(2.34) 


It  is  surprising  to  note  from  Eqs.  (2.33)  and  (2.34)  that 
the  leading  term  in  qx^/A  actually  adds  to  both  the  particle 
velocity  variance  and  its  mass  flux.  The  correlation  c'ui 
is  reduced  as  qx^/A  increases  but  the  leading  term  in  this 
reduction  appears  to  be  0(qx^/A)^  . 

Table  I gives  some  values  of  x^  for  different  size  par- 
ticles to  aid  in  determining  when  (qx^^/A)  <<  1 is  a valid 
approximation.  Since  A grows  proportional  to  z near  the 
surface  it  is  possible  to  estimate  the  height  at  which 
qxv/A  = 1 for  a given  typical  value  of  q(l  m/sec)  . At  large 
distances  above  this  height,  the  approximation  should  be  valid. 
Thus,  for  most  aerosols  smaller  than  that  associated  with  a fine 
mist  (a  = 50  pm)  the  approximation  is  valid  at  heights  greater 
than  a few  meters  above  the  surface. 

In  the  surface  layer  itself  we  must  consider  a more  com- 
plete analysis.  There  is  also  a problem  with  the  extremely 
small  particles.  If  the  stopping  distance  qx^  is  less  than 
the  Kolmogorov  microscale,  n , of  the  turbulence,  then  the 
particle  cannot  be  transferred  all  the  way  to  the  surface  by 
turbulence.  Instead  it  must  be  diffused  by  Brownian  motion 
across  the  surface  sublayer.  From  Tennekes  and  Lumley  (1972), 
n = (v^/e)^  and  from  our  model  e = q^/8A  , thus 

n = (8A)^/^(v/q)^/^  (2.35) 

The  minimum  q occurs  when  A ~ q at  the  edge  of  the  sub- 
layer, so  that 


q 


min 


2v/q 


(2.36) 


Table  I includes  values  of  q^x^/2v  which  show  that  for  tur- 
bulent air  with  q s:  i m/sec  typical  particles  with  diameters 
less  than  approximately  1 ym  must  involve  Brownian  diffusion 
through  a portion  of  the  surface  layer.  Table  I also  in- 
cludes values  of  D for  these  particles.  When  q^x.y./2v  > 1 , 
particles  will  be  carried  by  direct  inertial  impact  to  the 
surface.  For  particles  below  approximately  50  ym  the  approx- 
imation of  qxv/A  <<  1 will  describe  this  turbulent  diffusion 
down  to  the  level  we  would  ordinarily  apply  our  surface  boundary 
condition.  For  larger  particles  we  must  allow  qx^/A  to  be 
order  one  or  larger  in  the  surface  layer. 


The  time  scale  of  the  turbulence  in  the  surface  layer  is 
sufficiently  small  that  it  is  generally  close  to  equilibrium 
conditions  with  negligible  contributions  from  the  left-hand 
side  of  Eqs.  (2.23)  and  (2.28).  As  long  as  attention  is  restricted 
to  the  constant  shear  stress  region,  the  turbulent  diffusion 
terms  in  the  correlation  equation  also  contribute  little  to  the 
balance.  This  gives  validity  to  the  approximation  termed 
"Super -Equilibrium"  by  Donaldson,  1973.  The  Super-Equilibrium 
approximation  to  Eqs.  (2.20)  and  (2.23)  for  the  surface  layer  is 
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TABLE  I 

AEROSOL  TRANSPORT  PROPERTIES 

(Spherical  particles  with  * 1 gm/cm^  in  air  at  20°C,  1 ATM;  and  q “ 1 m/sec) 


Particle 

Diffusion 

Relaxation 

Terminal 

Height  at 

stopping  distance 

diameter 

Coefficient  time 

velocity 

which 

microscale 

2a 

D* 

q*T  /A=l 

q*T  /2V 

(Um) 

cm* /sec 

(sec) 

(cm/sec) 

(m) 

V 

0.1 

6.7  X 10"^ 

8.8  X 10~® 

8.6  X lO"^ 

2.9  X 10~® 

0.5 

6.3  X 10~^ 

1.0  X lO"® 

1.0  X lO"® 

3.3  X 10"® 

1.0 

2.8  X 10“^ 

3.6  X 10"® 

3.5  X 10"® 

0.12 

5 

4.9  X 10“® 

8.0  X 10~® 

7.8  X lO"® 

12.3  X lO"® 

2.7 

10 

2.4  X lO"® 

3.2  X 10"^ 

3.1  X lO"^ 

5 X lO"^ 

10.6 

50 

7.7  X 10"® 
-2 

7.58 

.012 

2.6  X 10® 

3 

100 

3.1  X 10 

-1 

30.3 

.048 

1.0  X 10 

500 

7.7  X 10 

757. 

1.18 

2.6  X 10® 

Tpriedlander,  1977. 
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-U.  w;  c ul  

A 3z  A^  9x. 


ir  - =1^  fir  - ^ 


c u.  -c  u. 

Ai  1 


- 0 (2,37) 


- Ifr  - <='  Hr  - 0 

J J 1 


(2.38) 


i 

To  keep  our  approximation  simple,  the  three  terms  ul  ^ 

_ — 3 p ’ g Q * ^ J 

(u^  -Uj)  and  (u^  -uj)  will  be  dropped 

since  they  all  go  to  zero  in  both  limits  of  -*■  0 and 

Ty  ->•  » . The  usual  model  of  the  last  term  in  Eq.  (2.38)  should 
be  modified  to  permit  it  to  approach  zero  as  x -^  «>  . An 
appropriate  choice  appears  to  tje  ^ 


7111= 

3Xj  a(1  + Kiqx^/A; 


(2.39) 


In  the  limit  of  x-^  ^ 0 , Eq.  (2.39)  gives  the  decorrelation 
term  as  inversely  proportional  to  the  turbulent  time  scale, 
but  when  the  slip  time  Xy  , is  larger  than  A/q  the  term 
becomes  inversely  proportional  to  x-y  . 

Rather  than  attempt  to  determine  sufficient  new  models  to 
obtain  u^  w^  and  u|w^  from  Eqs.  (2.21)  and  (2.22),  we  will 

approximate  both  as 

uj^  = upj  = u[wy  1 + K2(qx^/A)^]  (2.40) 

This  approximation  may  be  rationalized  by  noting  that 


t 

^i  J ^v 


(2.41) 


as  long  as  the  nonlinear  term  in  Eq.  (2.18)  is  ignored.  If  we 
also  approximate  u'  as  an  harmonic  oscillation,  then 
Eq.  (2.41)  leads  to 


^ = u'u^  = iP^y/  jl  + (qx^/A) 


(2.42) 


r 


IS 


Equations  (2.37)  through  (2.42)  lead  to 

f£  /[l  + K,(qT^/A)=  ] 


(2.43) 


and 


Aqc'u[p  ^ 

A l+KiqT^/A_  ’ 


u^w 


3c 


A U+K.qT„/AI  3 

These  two  equations  may  be  rearranged  to  give 


3u. 

1 


L I * , I 

z ~ ^ ^ 3z 


(2.44) 


c w 


± l£ 

Aq  3z 


c'wT  = - ^ 


^ w'w' 

Aq  3z 


l+KiqXq^/A 


l+(A4-Ki)qTv/A 
Ll+K2(qT.^,/A)2  J 


Equation  (2.46)  may  also  be  written  as 


(2.45) 


(2.46) 


c w^  = c w 


l-f-(A-f-Ki)qTq^/A 

Ll+Kiqx^/A 


(2.47) 


which,  for  qt  /A  <<  1 , agrees  to  0(qT  /A)  with  Eq.  (2.34). 
This  leaves  us  with  two  coefficients  K and  K,  to  determine 
turbulent  transport  of  particles  in  the' surface  layer  when 

qV^  ^ • 

In  sximmary,  our  proposed  approximation  for  turbulent  trans- 
port of  aerosols  is  based  on  writing  the  diffusion  equation 
for  c as 


a;:  -1 

l£  -h- h = s - ^ 

at  3x^  3x^ 


5 ^ 

p ax^ 


rrrr  fl+(^+K,)qx.r/^\ 

'"i  Vi+TC.qx^lt  ^ J_ 


(2.48) 


with  the  mass  flux  equation  modified  to 


ac^uT 

L 

3t 


+ u 


3c ' ul 
1 

■j  axj 


3c 


T— r ^'^i 


+ V 


a(qA 
c 3Xj 


3c^ 


I+KT^qr^^Tn*^  3Xj  “ ^ 


c uj 


i)  _ M ~ ~l 

3Xj  A (1-t-Kiqx^/A) 


(2.49) 
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1 


with  the  only  restriction  that  c <<  1 . Above  the  surface 
layer,  qx^/A  will  be  small  and  only  the  mean  drift  term  in 
Eq.  (2.48)  will  be  important.  While  in  the  surface  layer 
where  qt-^/A  can  become  of  0(1)  , Eqs . (2.48)  and  (2.49)  lead 
to  the  same  Super-Equilibrium  relations  as  just  given. 


In  combination  with  Eqs.  (2.48)  and  (2.49)  we  must  allow 
the  effective  height  zq  at  which  the  aerosol  concentration 
reaches  the  surface  valuS  to  be  a function  of  qx^  as  well  as 
a function  of  surface  characteristics.  In  fact  for  many  at- 
mospheric problems  it  appears  that  this  change  in  surface 
boundary  condition  will  be  the  major  influence  of  particle  size 
on  turbulent  transport  of  aerosols. 

Aerosol  surface  deposition  data  is  typically  presented  in 
terms  of  a deposition  velocity,  v^j  , defined  by  the  surface 
flux  divided  by  the  surface  concentration  (e.g.,  Chamberlain, 
1967) . The  actual  value  of  vj  depends  on  the  height  chosen 
to  measure  the  'surface'  concentration  due  to  the  possibility 
of  strong  gradients  near  the  surface.  Within  the  neutral  con- 
stant flux  layer,  Eqs.  (2.48)  and  (2.49)  reduce  to 


with 


and 


K -1^  + T gC  = VjC 
3z  V® 


d'-d 


K = 


L 


H-(A+Ki)qTv/A 


A = 6.5z 


(2.50) 


The  general  solution  to  Eq.  (2.50)  may  be  written  as 


VdCd 

c = + (const) exp 

V®  L 


When  the  boundary  condition  that 


(2.51) 


c = c^  at 


z = z. 


is  applied  to  Eq,  (2,51)  it  leads  to 


Vj  V , 

+ (1— ^ )exp 
Cd  T^g  x^g 


J 


-x^g  I K-Mz 


(2.52) 


As  X -►  0 , K -*■  0.22  qz  and 


c_ 

c 


T72q  ^^^d  ^ 


(2.53) 


the  classical  logarithmic  concentration  variation.  For  finite 
values  of  Xy  Eq . (2.52)  can  be  used  to  determine  z^  when  v, 
at  any  z^  is  knox^m. 
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3.  MODEL  CALCULATIONS 

The  general  A.R.A.P.  atmospheric  boundary  layer  model,  which 
includes  the  fluid  mechanics  of  mean  flow  dynamics  and  turbulence ; 
transport,  evaporation,  and  condensation  of  water;  and  atmospheric 
radiation,  has  been  exercised  for  a variety  of  situations.  These 
exercises  have  been  selected  to  demonstrate  the  power  of  the 
model  in  predicting  a variety  of  archetypal  atmospheric  flow 
situations  including  fog  episodes,  diurnal  variation  of  stratus 
over  the  ocean,  the  general  structure  of  the  trade-wind  boundary 
layer,  and  an  evaporatively  cooled  gust  front.  These  exercises 
were  selected  with  particular  attention  to  corresponding  experi- 
mental measurements  of  similar  or  identical  flow  situations. 

3 . 1 Marine  Boundary  Layer  Calculations 

A major  study  of  the  interaction  of  turbulence  with  atmo- 
spheric radiation  in  a class  of  fog  and  stratus  cloud  episodes 
is  described  in  Appendix  A which  is  a reprint  of  work  to  be 
published  in  the  Journal  of  the  Atmospheric  Sciences . In  ad- 
dition to  the  illustration  of  fog  and  cloud  f ormatio^n , develop- 
ment and  dissipation  for  various  situations,  we  present  in 
Appendix  A a theory  of  the  moist  atmospheric  surface  layer 
which  provides  a quantitative  classification  of  foggy  and  fog- 
free  surface  layers  dominated  by  turbulence.  We  further  in- 
clude, where  available,  observations  of  fog  and  stratus  extent, 
cooling  rates  and  liquid  water  content  and  show  favorable  com- 
parison with  the  A.R.A.P.  model. 

A full  comparison  of  the  A.R.A.P.  predictions  of  the  dynamics 
and  evolution  of  the  boxandary  layer  would  require  measurements 
of  the  fluid  state  at  an  initial  state  and  at  some  subsequent 
time.  Unfortunately,  so  far,  we  have  been  unable  to  obtain  such 
sequential  measurements.  We  are  thus  only  able  to  provide  a 
verification  of  the  consistency  of  the  A.R.A.P.  model  in  the 
sense  that  the  turbulence  co-existing  with  the  measured  mean 
profiles  at  one  time  is  well  represented  by  the  turbulence  cor- 
relations predicted  by  the  A.R.A.P.  model  for  the  same  mean 
profiles.  This  is  demonstrated  in  Appendix  B where  calculations 
for  a trade-wind  boundary  layer  which  has  received  some  experi- 
mental study  have  been  carried  out.  These  calculations  based 
on  the  data  summary  and  analysis  of  Lemone  and  Pennell  (1974) 
and  more  recent  GATE  data  also  provide  some  limited  comparison 
with  alternative  turbulence  models  based  on  direct  three-dimen- 
sional calculation  of  the  turbulence  with  modeling  of  the  subgrid 
turbulent  correlations  (Deardorff,  1972;  Sommeria,  1976). 

As  is  well  known,  direct  solution  of  an  ensemble  of  in- 
stantaneous Navier-Stokes  equations  for  the  turbulent  motions 
of  the  atmosphere  is  not  a feasible  or  practical  computation. 

A class  of  turbulent  models  which  attempt  direct  calculation  of 
the  large-scale  turbulence  but  which  do  not  compute  the  small 
scales  (of  the  order  of  the  grid  size  and  smaller)  are  known 
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as  subgrid-scale  turbulence  models  (Deardorff,  1972).  In 
these  computations  the  critical  modeling  centers  on  the  sub- 
grid turbulent  correlations  which  are  structured  within  the 
grid  length  and  time  scales  to  lie  in  the  so-called  "inertial” 
range  of  the  energy  spectrum  of  the  motion  being  computed.  Such 
calculations  must  of  necessity  be  fully  three-dimensional  and 
time-dependent  - and  hence  computationally  voluminous  - even 
when  the  mean  motion  is  principally  one-  or  two-dimensional  or 
steady . 

The  essential  predictions  of  both  calculations  under  these 
conditions  are  the  dynamics  of  evolution  to  and  the  final  slowly 
varying  state  over  a period  of  several  hours  following  initial- 
ization. Some  arbitrariness  necessarily  underlies  the  com- 
parison since  the  predicted  quantities  are  not  conceptually 
identical.  The  A.R.A.P.  second-order  closure  model  predicts 
the  ensemble -average  of  the  system.  The  calculation  of  Sommeria 
predicts  the  temperal  average  of  the  horizontal  spatial  average 
of  the  system.  Specifically,  for  atmospheric  scale  motions 
of  the  order  of  30  minutes,  Sommeria  calculates  at  5 s time 
steps,  averages  the  variables  over  the  horizontal  spatial  co- 
ordinate and  then  averages  these  quantities  over  4 minute  periods . 
These  4 minute  temporal  averages  are  probably  the  quantities 
most  closely  related  to  the  A.R.A.P.  ensemble  average. 

The  basic  features  of  this  comparison  are  as  follows. 

Below  the  clouds  and  in  the  mixed  layer  the  turbulent  momentum 
flux,  vertical  velocity  variance,  and  virtual  heat  flux  compare 
very  favorably  as  shown  in  Figs.  4 to  7 in  Appendix  B.  There 
is,  however,  a sharp  difference  in  the  average  density  of 
liquid  water  in  the  clouds  predicted  by  the  two  models  which 
in  turn  strongly  affects  the  turbulence  structure  in  the  cloud 
primarily  through  the  liquid  water  droplet  I.R.  radiation  inter- 
action. The  average  liquid  water  content  is  exceedingly  low 
in  the  subgrid  turbulence  calculation  (10“5  _ 10“O)  and  perhaps 
too  high  (for  the  scattered  cumulus  case)  in  the  A.R.A.P.  cal- 
culation ( 10"^).  Correspondingly,  there  is  about  a factor  of 
three  to  five  difference  in  the  turbulent  humidity  fluxes . As 
a result  of  the  higher  levels  of  water  in  the  A.R.A.P.  simulated 
clouds,  turbulence  production  due  to  cloud-top  cooling  is  en- 
hanced over  that  in  the  subgrid  turbulence  model  and  there  are 
relative  maxima  of  turbulence  in  the  clouds  as  well  as  in  the 
mixed  layer  below  the  clouds.  The  resolution  of  these  dif- 
ferences between  the  two  models  remains  to  be  given.  We  are  con- 
vinced that  the  simple  A.R.A.P.  averaging  algorithm  for  the 
liquid  water  content  considering  the  fluctuations  about  the 
saturated  state  gives  results  which  agree  within  approximately 
10  percent  with  those  for  the  more  complex  Gaussian  moment 
method  used  in  the  subgrid  turbulence  models  (Sommeria  and 
Deardorff,  1976)  . We  have  also  sha/7n  good  agreement  between 
the  A.R.A.P.  liquid  water  predictions  and  measurements  of  the 
liquid  water  in  fog  and  stratus  events  (see  Appendix  A) . The 
discrepancy  between  the  two  models  likely  stems  from  the  nature 
of  averaging  in  the  widely  disparate  states  of  scattered 
cumulus  and  continuous  stratus. 
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3 . 2 Gust  Front  Formed  by  Evaporative  Cooling 


As  a demonstration  calculation,  we  present  the  results  of 
a problem  for  which  both  nonhydrostatics  and  water  change  of 
phase  are  important.  Rain  falling  through  unsaturated  air  will 
cool  the  air  as  evaporation  takes  place.  This  air  with  its 
negative  buoyancy  can  then  form  a downdraft  which,  when  it  im- 
pinges upon  the  ground  spreads  out  to  form  the  local  gust  front 
often  associated  with  severe  thunderstorms.  We  have  made  cal- 
culations of  this  phenomena  for  NASA  (Teske  and  Lewellen,  1977) 
to  determine  the  windshear  and  turbulence  distributions  present 
because  of  their  potential  hazard  to  aircraft  operations.  In 
those  calculations  the  influence  of  evaporation  was  simulated 
by  introducing  a cold  downdraft  of  some  specified  temperature 
defect  at  a given  height  Zmax  with  vertical  velocity  . 

The  downdraft  is  then  permitted  to  fall  and  interact  with  the 
ground  to  form  the  outward  spreading  gust  front. 

For  the  present  calculation,  we  introduce  a downward  jet 
at  1.6  km  height,  as  shown  in  Fig.  3.1,  with  no  virtual  poten- 
tial temperature  defect  with  respect  to  its  surroundings.  We 
will  show  the  results  of  two  computations.  For  case  A,  the 
jet  is  supersaturated  with  a relatively  large  liquid  water 
content,  and  for  case  B the  jet  is  dry.  By  comparing  the 
results  of  the  two  cases,  it  is  possible  to  clearly  see  the  role 
of  evaporation  in  accelerating  the  resulting  gust  front.  It 
should  also  be  noted  that,  whereas  Teske  and  Lewellen  (1977) 
deals  with  an  axisymmetric  downdraft,  the  present  calculations 
are  for  t:TO- dimensional  downdrafts.  Other  conditions  such  as 
height  of  the  simulation  domain  and  the  boundary  conditions  on 
vorticity  and  velocity  there,  boundary  conditions  at  the  outer 
edge  of  the  domain,  and  surface  roughness  are  all  the  same. 
Initial  conditions  are  also  prescribed  in  a similar  fashion 
by  taking  the  boundary  values  of  vorticity  and  himidity  at 
Zmax  extrapolating  them  linearly  to  their  surface  values 

of  Hoj  “ 0 and  H = (0.36)Hs  . The  turbulence  is  initialized 
as  isotropic  with  qmax(”3  m/sec)  sufficiently  small  that  it 
can  be  dominated  by  turbulent  production. 


The  resulting  temperature  pattern,  as  it  evolves  in  time, 
is  shown  in  Fig.  3.2.  Contours  of  constant  virtual  potential 
temperature  are  shown  at  three  times,  t = 200,  500,  750  secs, 
for  case  A.  At  t = 0 , 6-^-0  everywhere  as  it  remains  in 

case  B.  The  maximum  virtual  potential  temperature  defect  is 
5.6°,  and  occurs  3C0  seconds  after  the  start  of  the  run.  The 
absolute  temperature  of  the  air  actually  rises  as  the  air  moves 
down  into  regions  of  higher  pressure  by  compression.  But  in 
order  to  stay  saturated  at  the  higher  absolute  temperature,  it 
evaporates  some  of  the  liquid  water  present  and  is  cooler  than 
it  would  be  if  it  descended  along  a dry  adiabatic  curve.  At 
the  later  times,  the  -1  contour  line  shows  that  the  same  shape 
front  as  evidenced  in  Fig.  3.2  of  Teske  and  Lewellen,  1977 
(reproduced  here  as  Fig.  3.3)  has  developed. 


o 

o 

o 

O 

o 

o 

o 

in 

o 

in 

o 

in 

U) 

CM 

o 

in 

CM 

^N  • Z 


w 


Figure  3.2b  Contours  of  constant  virtual  potential  temperature  for  case  A with  evaporation 
at  t = 500  secs  after  initialization.  (See  Fig.  3.2a  for  contour  notation.) 
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Figure  3.2c  Contours  of  constant  virtual  potential  temperature  for  case  A with  evaporation 
at  t = 750  secs.  (See  Fig.  3.2a  for  contour  notation.) 
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The  mean  velocity  contours  at  selected  times  are  shown  in 
Figs.  3.4  through  3.7  for  both  case  A and  B.  Case  B represents 
a stagnating  jet  without  any  influence  of  buoyancy.  Whereas 
the  vertical  velocity  in  the  center  of  the  downdraft  accelerates 
in  case  A with  evaporation,  shown  in  Fig.  3.4,  it  only  decelerates 
from  its  top  boundary  value  in  case  B,  shown  in  Fig.  3.5.  The 
difference  in  horizontal  velocities  as  shown  in  Figs . 3.6  and 
3.7  are  even  stronger.  Evaporative  cooling  allows  the  velocity 
for  case  A to  reach  a value  which  is  approximately  2 . 5 times 
larger  than  that  for  neutral  downdraft.  If  the  downdraft 
velocity  at  the  top  of  the  domain  were  smaller  or  the  height 
of  the  domain  larger,  this  difference  between  case  A and  B 
would  be  further  enhanced. 

The  total  velocity  variance  is  shown  in  Figs.  3.8  and  3.9 
for  case  A and  B respectively.  Due  to  the  combination  of  pro- 
duction by  higher  mean  velocity  gradients  and  by  buoyant  pro- 
duction, is  an  order  of  magnitude  larger  with  evaporation 
than  it  is  without. 

Contours  of  constant  values  of  liquid  water  content  are 
shown  in  Fig.  3.10  for  case  A.  The  smallest  value  contour  shown 
would  represent  the  edge  of  the  visible  cloud.  Such  a condensa- 
tion cloud  appears  different  from  that  marked  by  a tracer, 
such  as  dust  within  the  gust,  since  some  condensation  takes 
place  in  the  warm  air  which  is  lifted  by  the  advancing  gust 
front . 

One  feature  of  our  evaporatively  cooled  downdraft  which 
is  not  necessarily  observed  in  practice  is  the  cloud  reaching 
all  the  way  to  the  surface.  Since  our  model  does  not  currently 
permit  any  rain  out  of  liquid  water  content , it  remains 
available  for  evaporation  as  the  air  is  heated  by  descent. 

In  the  real  world  the  liquid  water  contained  in  relatively 
large  rain  drops  is  not  available  for  evaporation.  The  water 
drop  falls  through  the  air  before  it  has  time  to  completely 
evaporate.  Thus  when  precipitation  is  included  it  is  possible 
for  H to  exceed  Hs  all  the  way  to  the  surface  without  the 
cloud  reaching  the  surface.  Although  this  feature  may  have  a 
strong  influence  on  the  shape  of  our  simulated  "cloud,"  it 
should  have  a relatively  weak  influence  on  the  gust  front 
dynamics . 

Figure  3.11  shows  a picture  of  the  leading  edge  of  a gust 
front  as  marked  by  dust  (Turner,  1973).  When  the  simulated 
front  shown  in  Fig.  3.3  is  visualized  in  terms  of  temperature 
defect  intensity,  with  the  vertical  and  horizontal  dimensions 
shown  to  the  same  scale,  then  the  leading  edge  appears  as  shown 
in  Fig.  3.12.  The  strong  qualitative  similarity  between  Figs. 

3.11  and  2.12  is  readily  apparent.  The  quantitative  predic- 
tions also  appear  consistent  with  available  observations . 
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Figure  3.4b  Contours  of  constant  mean  vertical  velocity  for  case  A with  evaporation 
at  t = 500  secs  after  initialization. 
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Figure  3.5c  Contours  of  constant  mean  vertical  velocity  for  case  B without  evaporation 
at  t = 750  secs  after  initialization. 
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Figure  3.7a  Contours  of  constant  mena  horizontal  velocity  for  case  B without  evaporation 
at  t = 200  secs  after  initialization. 
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Contours  of  constant  mean  horizontal  velocity  for  case  B without  evaporation 
at  t - 500  secs  after  initializati nn  ^ 


Figure  3.7c  Contours  of  constant  mean  horizontal  velocity  for  case  B without  evaporation 
at  t = 750  secs  after  initialization. 
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at  t = 200  secs  after  initialization. 
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Figure  3.8b  Contours  of  total  velocity  variance  for  case  A with  evaporation 
at  t = 500  secs  after  initialization. 


at  t = 750  secs  after  initialization. 
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water  content  contours  for  case  A with  evaporation 
200  secs  after  initiation. 


Figure  3.10c  Liquid 


Figure  3.11  The  leading  edge  of  a gust:  front;  as  marked  by  dust  from  Turner,  1973. 


Figure  3.12  Model  simulation  of  the  leading  edge  of  a gust  front  when  visualized  in 
terms  of  temperature  defect  intensity. 


52 


4.  CONCLUSIONS  AND  RECOMMENDATIONS 

We  believe  that  substantial  progress  has  been  made  during 
the  last  year  toward  the  goal  of  providing  the  Navy  with  the 
capability  of  making  detailed  predictions  of  the  low-level 
distributions  of  pertinent  atmospheric  variables.  Although 
verification  tests  have  not  been  adequate  to  provide  detailed 
assessment  of  the  accuracy  of  our  model,  the  one-dimensional 
version  of  our  model  should  be  capable  of  making  useful  predic- 
tions over  regions  with  strong  spatial  homogeneity.  A copy 
of  this  program  has  been  delivered  to  the  Naval  Environmental 
Prediction  Research  Facility. 

During  the  current  year,  we  have  significantly  extended 
the  two-dimensional,  unsteady  model.  However,  additional  work 
is  required  to  make  these  features  completely  operational  when 
all  the  other  features  of  the  program,  such  as  density  gradients, 
condensation,  terrain  variation  and  thermal  radiation  are  being 
computed  simultaneously.  This  consolidation  of  the  model's 
capability  will  provide  a unique  research  tool  for  simulating 
much  of  the  physics  which  occurs  in  the  atmospheric  boundary 
layer  in  the  vicinity  of  coast  lines. 

As  demonstrated  in  Appendix  A,  we  have  been  able  to  use 
our  one-dimensional,  unsteady  model  to  simulate,  quite  realis- 
tically, many  of  the  features  observed  in  the  evolution  of 
coastal  fog  and  low-level  stratus.  More  features  should  be 
amenable  to  simulation  with  the  extended  two-dimensional  version 
of  our  model.  Ideally,  by  calculating  model  predictions  for 
situations  that  have  been  carefully  observed,  it  should  be  pos- 
sible to  both  make  model  verification  comparisons  and  suggest 
causal  relationships  between  different  physical  features  to  aid 
in  parameterization  for  larger  scale  models.  We  will  continue 
to  search  for  available  data  to  be  used  for  our  model  comparisons. 

Under  certain  conditions , the  large  eddies  in  the  marine 
boundary  layer  become  highly  organized.  Such  boundary- layer 
roll  vortices  have  been  the  subject  of  several  experimental 
measurement  programs  (e.g.,  LeMone,  1976).  By  fixing  our  com- 
putational coordinate  system  with  respect  to  the  axis  of  an 
individual  roll  vortex,  we  should  be  able  to  compute  the  large 
eddy  velocities  as  a part  of  the  ensemble-averaged  mean  velocity. 
The  numerical  results  may  then  be  compared  with  the  aircraft 
observations  published  by  LeMone.  In  addition  to  being  a basic 
flow  structure  of  interest  in  its  own  right,  such  a comparison 
between  model  results  and  field  observations  would  provide  a 
further  verification  test  of  the  model.  The  only  model  develop- 
ment required  for  this  calculation  is  the  application  of  periodic 
boundary  conditions  to  the  two-dimensional,  unsteady  program. 

We  propose  to  continue  our  investigation  of  the  structure 
and  movements  of  sharp  horizontal  temperature  gradients,  i.e., 
fronts,  in  marine  boundary  layer.  By  combining  an  analytical 
investigation  of  the  relationships  between  the  basic  parameters 
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in  a highly  simplified  model  of  the  front  as  an  internal  wave, 
with  detailed  calculations  using  our  sophisticated  two-dimensional 
numerical  model,  we  should  be  able  to  increase  our  understanding 
of  this  phenomenon.  We  will  look  for  frontal  situations  with 
characteristics  which  may  be  anticipated  to  have  an  impact  on 
naval  operations. 

Our  program's  usefulness  can  be  increased  by  decreasing  the 
computational  time  required.  We  have  taken  several  steps  in 
this  direction  over  the  last  year  by  changing  the  structure  of 
the  numerical  convective  operator,  by  adding  a quasi-coupling 
between  some  of  the  equations,  and  by  structuring  the  surface 
boundary  conditions  to  permit  larger  grid  spacing.  We  expect 
to  continue  this  process  of  improving  the  program's  efficiency. 

We  are  particularly  interested  in  achieving  a sufficiently 
significant  decrease  in  computational  requirements  to  make  the 
extension  to  a three-dimensional  model  attractive. 

The  current  analysis  of  aerosol  transport  lays  the  founda- 
tion for  incorporating  more  cloud  physics  into  our  model.  We 
would  like  to  pursue  this  far  enough  to  at  least  include  a 
simple  model  of  precipitation  in  our  program.  This  will  re- 
quire carrying  additional  information  about  particle  size  and 
growth. 
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THE  INTERACTION  BETV/EEN  TURBULENT  AND 
RADIATIVE  TRANSPORT 
IN  THE  DEVELOPMENT  OF  FOG 
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ABSTRACT 

A turbulent-radiative  description  of  the  low  level  atmo- 
sphere incorporating  stratus  cloud  and  fog  is  developed  based 
on  the  methodology  of  second-order  closure.  Existence  criteria 
for  surface  layer  fog  are  developed  on  the  basis  of  the  tur- 
bulent description.  Illustrative  application  is  made  to  the 
formation  and  structure  of  several  fog  and  cloud  episodes  in- 
cluding advectlve-radiative  fog,  subsidence  capped  stratus 
over  the  ocean,  and  surface  fog  resulting  from  the  nocturnal 
lowering  of  stratus. 

*To  be  published  in  the  Journal  of  the  Atmospheric  Sciences, 

Feb.  1978. 
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1.  Introduction 

In  this  study  we  examine  the  mutual  interaction  of 
turbulence,  radiation,  and  water  condensation  in  the  forma- 
tion, structure,  and  evolution  of  the  atmospheric  boundary 
layer  with  fog  or  clouds.  We  formulate  a turbulent-radiative 
theory  of  the  atmospheric  boundary  layer  based  upon  a second- 
order  closure  description  of  turbulence  which  includes  the 
condensation,  evaporation,  and  transport  of  v:ater  both  as 
vapor  and  liquid.  The  methodology  of  second-order  closure 
consists  of  extending  the  variables  which  describe  the  fluid 
beyond  the  mean  velocity,  humidity,  and  temperature  fields 
to  include  the  second-order  turbulent  correlations  of  these 
variables.  This  extended  set  of  correlations  is  governed  by 
an  extended  set  of  conservation  laws  obtained  by  taking 
exact  second-order  mioments  of  the  instantaneous  Navier-Stokes 
system..  Closure  is  obtained  by  m.odeling  the  unknovm  correla- 
tions which  appear  in  the  mom.ent  equations  in  term.s  of  the 
describing  field  of  second-order  correlations  (Donaldson,  1973; 
Mellor  and  Yamada,  1974;  Lumley  and  Khajeh-l^ouri,  1973; 

Wyngaard,  et  al.,  1974). 

The  second-order  closure  model  we  apply  has  previously 
been  used  to  obtain  accurate  prediction  of  the  ?^onin-Obukhov 
similarity  functions  for  the  surface  layer  (Lev/ellen  and 
Teske,  1973)  as  well  as  a detailed  representation  of  the 
diurnal  variation  of  the  atmospheric  boundary  layer  over  land 
(Lev;ellen,  Teske,  and  Donaldson,  1974),  and  the  sea-breeze-land  bree 
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cycle  in  the  coastal  planetary  boundary  layer  (Lewellen  and 
Teske,  1976).  A similar  second-order  closure  model  has  been 
applied  to  the  moist  boundary  layer  (but  neglecting  the  presence 
of  condensed  water  and  radiative  transport)  by  Burk,  1977- 
From  such  a turbulent-radiative  theory,  we  obtain  a 
description  of  the  dynamics  of  a cloudy  or  foggy  boundary 
layer.  At  the  scale  of  the  cloud  Itself  we  describe  its 
internal  dynamical  and  thermal  structure  including  the 
distributions  of  turbulent  heat  and  water  fluxes  and  radiative 
heating  and  cooling.  Of  particular  interest  at  this  scale 
is  the  detailed  behavior  of  cloud  top  mixing  in  which  radiative 
cooling  and  turbulence  production  are  strongly  linked.  We 
show  that  because  of  the  unique  radiative  properties  of  the 
water  molecule,  stratus  cloud  and  thick  fog  layer  tops  are 
not  "burned  away"  by  the  sun's  direct  radiation;  rather  the 
tops  of  such  clouds  and  fog  are  evaporated  away  by  turbulent 
transfer  of  solar  energy  captured  in  the  cloud  or  fog 
Interior, 

At  the  scale  of  the  overall  boundary  layer,  v;e  describe 
the  dynamics  of  formation  and  dissipation  of  cloud  and  fog 
events  subject  to  diurnal  or  mesoscale  boundary  and  initial 
conditions.  In  particular,  we  formulate  in  some  detail  the 
conditions  for  existence  of  surface  fog  in  a turbulent  surface 
layer  as  v;ell  as  the  conditions  for  diurnal  lowering  of 
stratus  to  the  surface. 

At  present  we  carry  out  this  program  of  description  and 
Illustration  v:ith  a one-dimensional  unsteady  and  two-dimensional 
steady  state  version  of  turbulent-radiative  theory. 

L ...  


V'e  thus 
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restrict  attention  to  layer  fog  and  stratus  cloud.  We  expect 
such  a turbulent-radiative  theory  based  upon  second-order 
closure  to  be  Instructive  in  the  development  and  evaluation  of 
simple  heat  and  moisture  exchange  models  betv;een  the  surface 
and  the  general  circulation  (Lilly,  1968;  Schubert,  1976; 

Randall,  1976).  Such  a theory  should  also  be  complementary 
to  alternative  approaches  to  detailed  models  of  atmospheric 
turbulence  employing  sub-grid  turbulence  modeling  (Sommerla  , 

1976). 

In  Parts  2 and  3 we  formulate  the  theoretical  model. 

I 

[ We  then  illustrate  its  descriptive  and  predictive  power  in 

[ application  to  surface  fog  existence  criteria  and  the  struc- 

^ ture  of  fog  in  the  turbulent  surface  layer  (Part  h) , advective- 

radiatlve  fog  (Part  5),  stratus  cloud  in  the  boundary  layer 
(Part  6),  and  diurnal  lowering  of  stratus  to  the  surface  as  fog 
(Part  7) . 

i 

I 


I 


A-5 


2.  Turbulent  description  of  the  atmospheric  boundary  layer 
The  moist  planetary  boundary  layer  can  be  described  by 
its  instantaneous  velocity  field  u^  , potential  temperature 
9 , total  water  mixing  ratio  q , and  pressure  p . These 

5 

variables  are  quite  accurately  governed  by  the  Instantaneous 
Navier-Stokes  equations  in  the  Boussinesq  approximation  in 
a coordinate  system  x^  = Xj^(x,y,z)  rotating  with  angular 
velocity  in  time  t : 


(la) 


^^i  ^ 1 ^ 

Dt  p 8x^ 


ijk^j'^k 


(lb) 


(Ic) 


reference  temperature  and  density  of  the  undisturbed  atm.osphere 
about  which  the  Boussinesq  approximation  is  made.  The 
source  term  for  radiant  energy  is  S and  m.ay  be  expressed 
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Our  description  is  restricted  to  ncnprecipitating  boundary 
layers  so  that  there  is  no  source  term  for  moisture  in  Eq. 

(Id).  Consistent  with  this  approximation,  vre  will  restrict 
this  analysis  to  equilibrium  condensation  so  that  the  vapor 
mixing  ratio  and  the  liquid  mixing  ratio  q^^  are  given  by 

‘^v  " ^ 

= (q-qg)^(q-qs)  (3b) 

where  ff( • ) is  the  Heaviside  function  and  q_(p,T)  is  the 
saturation  mixing  ratio. 

The  conserved  potential  temperature  0 is  related  to 
the  absolute  temperature  T as 

e.,  = T - + fz  + ^ q (H) 

s r ^p 

where  F = S/Cp  is  the  adiabatic  lapse  rate,  L is  the 
latent  heat  of  vaporization,  and  z is  the  coordinate  along 
which  the  gravitational  vector  g^  lies.  The  virtual  poten- 
tial temperature  9^  is  the  active  temperature  in  the 
buoyancy  force  defined  as 

= T - T^  + Fz  + (0.61q^  - q^)T^  (5) 

The  Dotential  temoeratures  9 and  9,  are  necessarily  eaual 

* “ S V - 

in  a vrater-free  atmosphere  but  not  in  a water-laden  boundary 
layer.  If  all  X'^ater  is  in  vapor  form  so  that  q_^  = q 


then 
or  9g 


1 


! 


[ 

i 

i 

i 


the  energy  equations  expressed  either  ir 


erns  of  9 
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The  parameters  y , u , a , S are 


y = 1 + 1.61  a = 1 + 1.61  ^ 


'll'.". 
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3£nq^ 

a£nT 


(9) 


If  the  saturation  function  is  given  in  the  Clape^/ron  representa- 
tion 

(10) 


(p“)  _R^  (t„  - t)_ 


where  , p„  , T„ 

then  is  given  by 


are  reference  values  at  some  state 


Bt  = R^T 


(11) 


Equations  (la,b,d)  and  (4)  form  the  basis  of  the  turbulent 
description.  We  now  consider  their  turbulent  ensemble  average. 
Denote  the  ensemble  average  with  an  cverbar  and  let  a prime 
denote  the  deviation  from  the  average.  The  mean  fluid  equations 
for  u^.  , 5^  , q , are  then 

3u< 


= 0 


* i 

Dt  p 3x^ 


^ t;  ®v  - 


D0 


V 


Dt 


(u'. 0')  + < s + < rw 
0x.  ^ J v'  s g 


Dq 

Dt 


(12a) 

(12b) 

(12c) 

(12d) 


In  Eqs.  (12)  it  is  understood  that  the  substantial  derivative 

is  based  uoon  the  mean  velocity  u,  . In  the  ensemble 
Dt 

averaging,  we  neglect  fluctuations  in  the  coefficients 


T 
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\'Je  nov;  augment  the  fundamental  set  of  mean  flov;  variables 
to  include  the  following  second-order  correlations; 


uTuT 
i j 


(Reynolds  stresses) 
(heat  flux) 


(moisture  flux) 


o'e ' 
H V 


( temperature-humidity 

correlations  and 
variances ) 


These  second-order  correlations  are  governed  by  an  exact  set 
of  moment  equations  generated  from  the  instantaneous  Navier- 
Stokes  system  (la-d)  (Donaldson,  1973;  Lewellen,  1977). 
Unknovm  second-order  correlations  and  third-order  correlations 
which  appear  in  the  exact  moment  equations  are  modeled  in 
terms  of  second-order  correlations  leading  to  the  follov/ing 
closed  system  (Lev;ellen  and  Teske,  1975) 


D -7— r 

Dt  ^i^J 


£,•  g. 

"^l^k  ~ '^J’^k  9177  ^ ^ ^’'®v 


- 2 = S,“J  - “k“i 

^ / 9uTuT  \ 


qjQV3) 


ij  12A 


(13a) 


Dt  V 


3x. 


1 V 3x. 


8ui9, 


- 2-iJk^J  W " °*3 


- 0.75  ^ up;  + <^T  u^w’ 


D 


3u. 


-ujuj  - u:^'  ^ ^ q'e:. 


13^— 13^^  '' 

8 


9u 


" ^k*?  ' 3x^ 


a A 


j 


3x , 


- 0.75  ^ n[q' 


D 


39. 


Dt  ^ 


= -u 


TTf 


3o 


\3x. 


- < r 6 


ij 


+ 0,3 


ax. 


9<?  '0 ' \ 

-33r^  - °-'5  ipp: 


j 


5_  flt2  _ 

Dt  ®v 


ae. 


°-3  air  Ip 


- 0.45  ^ + 2k  r 

A V g V 


D — 
Dt 


“■3  jIj  (=A 


- 0.45  ^ ^ '2 


(i; 


(i; 


(i; 


(i: 


(i; 


In  the  development  of  Eo.s . (13)  correlations  of  the  type 
3 ' u|  , s ' i and  S'q'  appear.  These  radiative  correla- 
tions arise  mainly  because  of  fluctuations  in  the  radiation 


absorption  coefficien 


for  v;ater  droplets  vhich  in 


turn  arise 
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because  of  the  fluctuations  in  . We  observe  that  these 

correlations  may  modify  the  coefficients  of  the  turbulent 
decay  terms  in  Eqs.  (13b, d,  and  e);  but  their  structure  is 
not  likely  to  be  modified.  In  the  present  calculation, 
we  do  not  modify  the  coefficients  of  the  model  in  Eqs.  (13) 
to  account  for  possible  differences  in  radiative  and  molecular 
diffusion  in  small  eddy  destruction,  since  such  effects 
require  separate  empirical  conf irm.ation  and  we  do  not  expect 
the  modifications  to  be  large. 

The  variable  q appearing  in  all  correlation  equations 
is  the  rms  turbulence  velocity.  We  note  that  each  of  the  cor- 
relation Eqs.  (13)  is  a balance  equation  involving  production/ 
destruction,  diffusion,  and  turbulent  decay.  Three  distinct 
production/destruction  mechanisms  are  operative  in  the  planetary 
atmosphere  corresponding  to  the  first  terms  on  the  right-hand 
side  of  the  correlation  equations.  These  three  mechanisms 
are  raea.n  variable  gradients,  buoyancy,  and  rotation  (Coriolis 
forces ) . 

Central  to  all  turbulent  decay  and  diffusion  processes  in 
the  correlation  equations  is  the  turbulent  scale  A . This 
scale  function  must  also  be  determined  to  complete  closure. 

The  form  of  the  scale  equation  is  similar  whet;hsr  one  begins 
with  the  tv;o  point  velocity  correlation  to  define  an  integral 
scale,  the  dissipation  function  c = q"/3A  , or  the  vorticity 
function  q/A  . In  all  cases  the  coefficients  of  the  individual 
terms  must  be  determined  empirically.  The  form  of  the  scale 
equation  used  in  the  present  v/crk  is  (Lewellen,  1977) 
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eA  = 0.35  ^ ^ + 0.75,+  0.3  ^ 

Q <3  «3 

J ^ ^ 


(•■%) 


The  correlations  ql  , q' q'  , and  q' q'  (v/hich  do  not 

S S V s 

appear  explicitly  in  the  moment  equations)  are  also  useful. 

T-l IT  I I ^ 1__  i 1 


For  small  fluctuations 


in  terms  of  the  basic  set  as 


<<  1 these  may  be  represented 


2 —2  Mr 


d'~  + 2 ^ y q'Q'  + y <?  ' 
2v  ^ 


^ ■ VsM  \ ‘I'K 

L r ^ 

=3x^3  ^ ^Tq-T  + y 


where 


^0  = 


-.61 


and  upper  and  lov;er  values  apply  to  nonsaturated  and  saturated 
cases,  respectively. 

V/e  now  turn  to  the  specification  of  the  ensemble  average 
liquid  water  content  which  is  given  by 

^ (19) 

This  average  is  not  equal  to  the  simple  form 

^ (20) 

This  is  because  the  liquid  water  content  q^  is  a physically 
positive  definite  quantity.  It  may  be  seen  that  the  contributions 


r 
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made  when  q+q'  > q +q'  cannot  be  balanced  by  contributions 

s s 

made  when  q + q'  < Q +<? ' • Thus,  the  form  Eq.  (20)  is  correct 

s s 

only  when  the  flow  is  sufficiently  saturated  that  no  local 

unsaturated  region  is  contained  within  an  eddy. 

The  exact  specification  of  the  average  liquid  water  ratio 

may  be  calculated  by  taking  the  moment  implied  by  Eq.  (19) 

over  the  Joint  probability  distribution  function  f(f?,<?_) 

s 

for  q and  q^  . Since  the  statistical  detail  of  second-order 
closure  is  limited  to  second-order  moments,  the  distribution 
function  fiq,q  ) is  not  knov/n  or  describable.  Only  the 
second  moments  of  f may  be  elements  of  the  theory.  The 
liquid  water  moment  may  be  carried  out  by  assuming  the  form 
of  ^(^jq-)  e.g.,  a joint  Gaussian  distribution  about  the 

mean  (Sommeria  and  Deardorff,  1977;  Mellor,  1977),  however 
it  is  possible  that  quite  arbitrary  artificiality  may  be  in- 
jected into  the  theory  at  significant  formal  expense  by  such 
assumed  distribution  functions.  Vie  adopt  an  equally  arbitrary 
but  much  simpler  procedure  for  the  determination  of  q^ 
which  requires  assumptions  about  the  nature  of  the  averages 
rather  than  the  nature  of  the  distribution  function. 

To  account  for  the  fact  that  eddies  may  contain  both 
saturated  and  unsaturated  fluid,  we  define  a transition  regime 
which  separates  completely  unsaturated  and  completely 
saturated  regimes.  This  occurs  when  the  majority  of  the 
fluctuations  in  q and  q give  q-q^  and  q-q„ 

o So 

opposite  senses.  This  transition  regime  m.ay  be  defined  as 
occurring  v/henever 
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(21) 


where  A is  a characteristic  measure  of  the  size  of  fluctua- 
tion which  may  contain  both  saturated  and  unsaturated  eddies: 


(22) 


We  then  approximate  the  average  liquid  water  content  as 


|c^-«3+a] 


A < U-<?3l 
A > U-Qgl 


(23) 


Equations  (23)  have  the  correct  limits  of  zero  and 
v/hen  equality  holds  in  Eq.  (21).  Further,  if  q'  and  q' 
are  each  symmetrically  distributed  about  their  means,  Eq.  (23) 
indicates  that  when  q = and  <?  ’ (or  q')  vanishes,  only 

S o 

negative  q'  (or  positive  q' ) will  lead  to  condensation, 
s 

The  coefficients  k , k also  reauire  modification 

s ’ g 

in  the  transition  regime.  We  structure  these  coefficients 
in  the  transition  regime  to  vary  linearly  between  their 
saturated  and  unsaturated  values  in  direct  proportion  to  the 
liquid  water  content.  Thus,  in  the  transition  regime  defined 
by  the  inequality  (21) 


■=S  ■ 1 * (?  - 

^ ^25) 


otherwise,  k and  < are  given  by  Eqs. 
s g 

V/e  nov/  remark  on  boundary  conditions 
system  consisting  of  Eqs.  (12)  and  (13). 
diticn  on  the  mean  velocity  at  the  surface 


(7). 

appropriate  to  the 
The  boundary  con- 
is  given  in  term.s 


T 
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of  the  effective  aerodynamic  roughness  . Over  v/ater, 

the  effective  roughness  height  is  given  by  Froude 

similarity  scaling  (Charnock,  1955).  k'e  use  the  Froude 
numiber  based  on  u*  determined  em.pirically  by  V/u  (1969) 


I 

i 

I 

i 

i 

t 


— it-  = 8 (26) 

/gi^ 

Surface  boundary  conditions  are  applied  at  a reference 
height  z = which  is  greater  than  z^  but  sufficiently  small 

that  two  conditions  are  satisfied.  First  the  characteristic 
time  of  the  dynamics  of  turbulence,  A/q  , is  small  enough 
so  that  the  turbulent  quantities  are  in  a local  equilibrium 
balance  between  production  and  decay  and,  second,  that 
falls  v;ell  within  the  constant  flux  region  where  the  turbulent 
flux  of  momentum,  heat,  or  species  remains  constant  with 
height.  The  boundary  conditions  applicable  at  z^  are  that 
the  gradients  of  mean  flow  variables  are  related  to  the 
surface  fluxes 


0 

= w'9' 

, Q:i  = y<'  Cl' 

o ’ 

vs 

V 

Q'  * 3 -w  ^ 

through  the  Monin-Cbukhov  similarity  functions  (Lewellen  and 
Teske,  1973): 
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and  ? is  the  height  z nondimensionalized  by  The  Monin- 
Obukhov  length.  The  similarity  functions  and  fg  have 

been  determined  empirically  (Businger,  et  al.,  1970)  and 
theoretically  (Lewellen  and  Teske,  1973;  Mellor,  1973)* 

Over  the  ocean  surface,  values  for  H and  9^  are 

more  likely  to  be  known  than  the  surface  fluxes  w'q' 

o 

w*9'  . In  this  case,  Eqs . (27)-(29)  are  analytically 

^ o 

Integrated  between  z^  and  and  the  surface  values 

®vo  ’ ^o 

At  the  upper  boundary  of  the  computational  domain  (which 
lies  outside  the  boundary  layer)  we  specify  the  lapse  rate, 
geostrophic  wind,  and  upper  air  v;ater  mixing  ratio.  All 
turbulence  quantities  are  made  to  vanish  at  the  upper  boundary 
for  stable  atmospheres.  Notice  that  the  lapse  rate  specifica- 
tion allows  the  boundary  temperature  to  be  computed  in  the 
course  of  the  dynamics  of  evolution  of  the  boundary  layer. 

This  is  particularly  important  because  a diurnal  radiative 
cooling/heating  in  general  exists  on  the  boundary  as  well  as 
in  the  interior  of  the  fluid. 

Sample  calculations  using  this  system  of  equations  and 
boundary  conditions  were  reported  earlier  (Lev;ellen  and 


I 
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T 


Teske,  1975)  without  including  the  influence  of  coupling 
to  thermal  radiation. 


T 
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3.  Radiant  heat  transfer  description 

We  separate  the  direct  solar  (short  wavelength)  from 
the  terrestial  (long  wavelength)  fields.  Each  of  these 
fields  Is  described  in  terms  of  the  one-dimensional  (z) 
upward  and  downward  frequency  and  angular  averaged  radiation 
Intensities  f"*"  , F”  . These  intensities  are  governed  by 
the  upward  and  downward  transport  equations  in  integral 
form  (Goody,  1964). 


F‘^(z) 


[F'‘'(z^)  - a(zQ)]T(z^,z)  + a(z) 
z 

- I t(z’  ,z)  ||-r  dz' 


(31) 


F“(z)  = [F“(h)  - a(h)]T(h,z)  + a(z) 
z 

- t(z',z)  ||-r  dz’  (32) 

h 

The  boundary  values  at  the  surface  z = z^  and  boundary 
layer  "top"  z = h are  F'*’(z^)  and  F~(h)  respectively. 

The  radiative  heat  flux  given  in  terms  of  F"*"  , F“ 

for  either  radiation  field  as 


(33) 


and  the  total  radiant  heat  flux  is  the  sum  of  the  direct  solar 
heat  flux  and  the  terrestial  heat  flux. 

The  absorption,  scattering,  and  emissive  properties  of 
the  atr.osphere  are  embedded  in  the  source  function  a(z)  and 
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the  transmission  function  tCz^jZ^)  between  two  levels  z^,  Z2  . 
We  neglect  scattering  in  the  present  analysis.  For  direct 
solar  radiation,  emission  is  neglected  and  a(z)  =0  . The 
boundary  intensity  F~(h)  is  set  equal  to  the  incident  solar 
flux  at  h while  F'*’(z^)  is  determined  by  the  surface 
reflectivity.  The  transmission  function  t is  modeled  after 

that  of  Manabe  and  Strickler  (1964)  for  water  vapor  and  CO2 . 

To  account  for  liquid  water,  we  include  a transmission 

T^^q(z^,Z2)  as  a factor  in  the  total  transmission  t of  the  form 

’^liq^^l’^2^  = exp  (-3a®m)  (34) 

where  m is  the  absorbing  mass 

Z2 

m = pCdz  (35) 

^1 

6 is  the  diffusivity  factor,  and  is  the  specific  absorp- 

tion coefficient  for  direct  solar  radiation  by  water  droplets. 
Feigel'son  (1964)  notes  that  solar  absorption  is  not  sensitive 
to  droplet  size  for  wavelengths  less  than  2u  . Integrating 
the  spectral  liquid  water  absorption  coefficient  tabulated  by 
Feigel'son  over  the  direct  solar  spectrum,  we  find 

= 16  cm^/gm. 


For  the  terrestial  field,  the  source  function  a is  set 


equal  to  the  integrated  Flanck  function.  For  boundary  condi- 
tions, we  set  F''’(z^)  equal  to  the  surface  Flanck  function 
and  F“(H)  equal  to  the  terrestial  downward  flux  emitted  by 
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T^iq  = exp(-6a’^m)  (36) 

T 

Here  a Is  the  averaged  liquid  water  absorption  coefficient 
for  terrestial  radiation  which  depends  strongly  upon  the 

2 

droplet  size  distribution  ranging  in  value  from  200-1700  cm  /gra 

(Feigel'son,  1970)  for  spectra  with  average  drop  size  radii 

ranging  from  ^^.5  to  7 U.  We  m.aintain  the  droplet  size  effect 

T 

on  radiative  absorption  as  a parameter  through  a . Detailed 

Inclusion  of  variation  with  droplet  size  requires  prediction  of  the 

droplet  number  density  within  a cloud  or  fog  as  well  as  the 

liquid  water  content.  Experimental  evidence  (Mack,  et  al., 

1974)  does  not  validate  the  simplifying  assumption  of  uniform 

droplet  density  in  space.  Their  observations  indicate  an 

order  of  magnitude  variation  in  droplet  number  density  through  the 

boundary  layer  while  mean  droplet  size  varies  by  about  30 

T 

percent.  A characteristic  value  of  a for  fog  and  stratus 

2 

(mean  drop  radius  approximately  6 p ) is  600  cm  /gm. 

The  significant  difference  between  the  absorption  co- 
s T 

efficients  a and  a is  noteworthy.  The  absorption  co- 
efficients may  be  represented  in  terms  of  absorption  lengths  as 


'rad 


(37) 


and  exhibited  in  Table  I for  various  (spatial)  average  liquid  water 
contents  • These  values  indicate  that  water  is  radiatively 

more  active  with  the  terrestial  spectrum,  then  the  direct  solar 
spectrum,  for  drops  in  the  4-10  p range  typical  of  fog  and 
cloud.  This  sharp  difference  in  radiative  em.ission/absorption 


i 

I 

I 


T 
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in  the  terrestial  and  solar  spectrum  will  influence  strongly 
the  structure  of  cloud/fog  cooling  and  the  mechanism  by 
which  direct  solar  radiation  evaporates  away  clouds  or  fog. 

We  discuss  these  implications  in  Parts  6,  7,  and  8. 

Boundary  conditions  for  terrestial  radiative  transport 
are  P~  = 0 at  the  upper  boundary  and  P''’  = Bo  (3q  Is  the 
black-body  flux)  at  the  surface.  Por  direct  solar  radiation, 
P~  is  specified  as  the  incident  downward  solar  flux  at  the 


upper  boundary  and  P"*"  = 

tF“  at  the  surface 

where  t is  the 

surface  reflectance.  The 

upper  boundary  for 

radiative  transport 

is  set  at 

the  atmospheric 

scale  height  and  the  average  liquid 

water  content  and  average 

water  vapor  content 

of  the  atmosphere 

contained 

within  the  radiative  boundary  and  the  fluid  upper 

boundary 

is  specified. 

Table  I. 

Radiative  absorption  depths  for  fogs  & 

stratus 

j,s,T  _ 

1 

rad  s 

,T^ 

a 

<ci^> 

Average 

Direct  Solar 

Terrestial 

Liquid  Water 

Radiation 

Radiation 

Contents 

rp 

<q^>  (g/kg) 

,(m) 

rad 

Tenuous 

Fog 

. 025 

20,000 

600 

.1 

5,500 

167 

.2 

2,700  i 

rrn 

A 

1,350 

41 

Dense 

.8 

700 

21 

Stratus 

Tops 

1.0 

500 

16 
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Fcg  in  turbulent  dominated  surface  layers 
VJe  now  take  up  the  first  of  several  illustrative  cases 
based  upon  the  turbulent-radiative  model  we  have  .just  developed. 
We  first  consider  the  case  in  which  turbulent  transport 
dominates  radiative  transport  in  a surface  fog.  This  represents 
a simple  application  of  the  turbulence  theory  and  also  allows 
a useful  extension  of  the  Monin-Obukhov  analysis  of  the  surface 
layer  to  Include  fog.  The  surface  layer  of  the  boundary 
layer  is  that  region  in  v/hich  the  turbulent  fluxes  of  v/ater, 
momentum,  and  heat  may  be  taken  as  constant . 

In  the  absence  of  radiation  and  for  fog  layers  v;ithin 

the  surface  layer,  the  turbulent  correlations  for  6 and  q 

s 

are  Identical  to  one  another  (although  dissimilar  to  the 

correlation  equations  for  u^).  The  distributions  of 

and  G -0  as  well  as  their  turbulent  correlations  must  there- 
s so 

fore  be  similar  to  one  another  and  a function  of  and 

z/z  only.  Surface  sim.ilarlty  leads  to  the  formation  of 
c 

the  following  theorem 
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as  essentially  noted  by  G.  I.  Taylor  (1917): 

All  fluid  states  in  a surface  layer  in  turbulent 
interaction  must  map  onto  a straight  line  in  the 
q-6g  plane  v;hich  passes  through  the  surface 
state  (c^Q,  ©sq)  and  has  a slope  determined  by 
the  surface  humidity  to  heat  flux  ratio  <<s/®s*  • 

Thus,  surface  layer  similarity  provides  the  basis  for  existence 

criteria  of  turbulent  dominated  surface  fogs.  In  the  <J-0_ 

s 

plane  of  Figure  1,  the  surface  state  is  denoted  by  "0"  and 

the  straight  line  OCDE  represents  the  distribution  of  q and 

6 in  the  surface  layer.  The  saturation  function  9„(p,T) 
s s 

may  be  represented  in  ^-0_  coordinates  and  is  also  shown 

s 

in  Figure  1.  [Over  the  range  of  the  surface  layer  (z  <<  T^F”^) 

we  neglect  the  variation  of  q with  pressure.]  The  tangent 

to  the  saturation  curve  at  the  surface  is  denoted  3^  and 

is  the  thermodynamic  derivative  with  respect  to  temperature 

in  q -0  coordinates: 
s s 


P 


The  slope  is  sim.ply  related  to  the  logarithmic  temperature 

derivative  Brp  as 
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For  saturated  conditions  at  the  surface,  any  fluid 
trajectory  will  be  a straight  ray  emanating  from  the  surface 
point  Qq  , . Such  a fluid  trajectory  entering  the 

region  below  the  line  with  slope  6^  will  be  a fog-free 
surface  trajectory.  This  region  is  determined  by  the 
saturation  conditions  which  determ.ine  6^  and  the 
turbulent  surface  layer  conditions  which  determine  j • 

Thus,  surface  fog  is  not  possible  for 


< 0 


e > 0 

s* 


s» 


ss 


> 6. 


(40) 


In  the  region  between  the  lines  v/ith  slopes  6^  and  , 

surface  fog  exists  adjacent  to  a cold  surface  and  at  a 
finite  depth  as  shown  by  the  typical  surface  trajectory 
0-C-D-E,  the  top  of  fog  being  at  D and  the  maxim.um  liquid 
content  at  C . Above  the  line  with  slope  and 

the  q axis,  cold  surface  fog  again  exists,  but  it  cannot 
be  bounded  within  the  surface  layer. 

In  the  region  of  9^^  > 0 , warm  surface  fog  exists.  Such 
a fog  may  be  bounded  within  the  surface  layer  with  a trajectory 
completely  analogous  to  0-C-D-E  for  cold  surface  fog  or  it  may 
be  unbounded  within  the  surface  layer. 
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V/e  observe  that  the  existence  or  nonexistence  of  surface 
fog  may  be  assessed  directly  in  terms  of  the  surface  flux 
ratio  qj./9  . and  that  the  critical  value  of  is 

equal  to  the  tangent  of  the  saturation  curve  at  the  surface: 
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nearly  adiabatic  air  over-running  a cold  water  surface  fogs 
easily.  On  the  other  hand,  for  a warr.  surface  fog  (0  > 0), 

the  flux  ratio  must  be  smaller  in  magnitude  than  the  critical 
value.  Hence,  cold  air  v/ith  high  relative  humidity  over- 
running a v/arm.  surface  v;ill  tend  to  fog.  Further,  for  fluid 
trajectories  with  ^s/0  ..  in  the  fog  regim^e  the  liquid  water 
distribution  is  determiined  by  the  difference  between  the 
straight  fluid  traj’ectory  OECD  and  the  saturation  curve.  In 
particular,  the  m.aximum  liquid  water  occurs  at  C and  the 
total  height  of  the  fog  layer  occurs  at  D in  Figure  1. 

An  illustrative  family  of  fog  distributions  as  given  by 
integrating  the  similarity  functions  from  Eqs . (28)  and  (29) 
is  shown  in  Figure  2.  of  particular  interest,  is  the  structure 
of  the  liquid  water  profiles  in  the  fog  layer  for  stable 
(C  > 0)  and  unstable  (?  < 0)  regir.es.  Stable  fog  layers 
have  miaximum  close  to  the  surface  and  drop  off  sharply 

in  liquid  water  content  up  to  their  m.aximum.  depth.  Unstable 
layers  are  well-m.ixed  and  do  not  shov;  a pronounced  maximumi 
except  for  very  thin  layers  <<  1 . Beth  families 

of  fogs  deepen  and  their  m.axim.um.  liquid  water  contents  increase 
as  their  flux  ratios  depart  m.ore  from  the  critical 

' S 7C 

value  B 

o 


Let  us  nov;  determ.ine  under  what  conditions  the  surface 
fog  layer  may  be  dominated  by  turbulent  rather  than  radiative 
transport.  The  surface  turbulent  heat  flux  will  be  determ.ined 
by  the  temperature  gradient  in  the  surface  layer;  but  this 
gradient  will  also  determ.ine  the  m.axim.um.  liquid  water  developed 
in  the  fog  layer  and  hence  the  radiative  emiission  from,  the 
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TURBULENT  DOMINATED  FOG 


Existence  diagram  for  turbulent  dominated  fog.  Slope  of  the  saturation  function 
in  Oq  space,  Bq,  determines  the  critical  q,.,/0s,v  surface  layer  distributions 

OCDE  which  will  be  foggy  or  fog-free.  For  a typical  foggy  layer,  fog  bank  height  i 
at  D and  maximum  liquid  water  content  at  C . 


DOMINATED  SURFACE  LAYER 


Illustrative  liquid  water  distribution  in  stable  (c  > 0)  and  unstable  (c  < 0)  surface 
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For  a radiatively  thin  fog  layer 
possible  to  estimate  the  radiative  to  turbulent  heat  flux 


for  the  fog  layer  as 


„rad 

a 


turb 

Q 


(43) 


For  the  nominal  condition  of  a 300°K  surface, 
and  a*^  “ 500  cm^/gm  , 


.00002 


. 6 £n(6/zn)  ... 

turb  2000  u. 
q 

a result  v;hich  depends  only  upon  the  strength  of  the  v:ind 

and  the  depth  of  the  fog.  We  see  that  it  takes  a relatively 

strong  v/lnd  to  make  negligible  compared  to  q^^^^ . 

_3 

For  6 = 10  m , = 10  m , a surf'ace  shear  velocity  of 

u*  = 0.6  m/sec  (corresponding  to  a neutral  wind  at  10  meters 
of  lA  m/sec)  or  greater  is  required  to  make  q^^'^  less  than 
10  percent  of  . We  conclude  that  radiation  plays  a 

significant  role  in  most  surface  fogs  of  depths  greater  than 
a few  m.eters. 
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5.  The  boundary  layer  with  advective-radiatlve  fog 

We  return  to  the  full  turbulent -radiative  description 
[Eqs.  (12),  (13),  (31),  (32)]  to  illustrate  the  more  general 
case  of  fog  formation  with  both  radiant  and  turbulent  trans- 
port. We  wish  to  utilize  that  description  to  predict  the 
formation  and  structure  of  a commonly  observed  class  of 
coastal  fogs  over  water  in  which  convection  and  radiation 
play  equally  important  roles.  Fogs  are  frequently  found  to 
occur  in  coastal  regions  subject  to  strong  upwelling  and  non- 
uniform  surface  temperatures  v/here  winds  bring  the  surface 
air  fromi  a region  of  cold  water  onto  a region  of  v/armer  water 
(warm,  surface  fog)  . Such  fogs  have  been  reported  and  documented 
by  Mack,  Pilie’,  and  Koclmtond  (1973),  and  Leipper  (1948).  As  an 
example  of  such  advective-radiative  fogs  predicted  by  the 
present  model  we  consider  the  following  archetypal  situation. 

A stable,  steady  boundary  layer  running  over  constant  temperature 
vfater  is  Incident  on  a region  of  linearly  increasing  surface 
temperature.  As  a result,  an  unstable  mixed  layer  will  be 
initiated  and  fog  may  form  in  the  nev/ly  developing  layer. 

The  incident  profiles  for  this  illustration  are  characterized 
by  an  = .00015  under  conditions  where 


The  full  profiles  are  determined  by  wind  velocity  of  2.8  m/sec 
and  an  absolute  temperature  difference  of  +1°C  at  z = 20  m 
elevation  above  the  sea  surface  of  T = 13°C  . Thus,  the 
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incident  flow  is  in  the  cold  surface  fog-free  regime  (Fig.  1). 
Direct  solar  radiation  was  suppressed  (simulating  nighttime 
conditions).  The  surface  humidity  v/as  maintained  at  its 
saturation  value  at  the  surface  temperature.  Surface  tempera- 
ture linearly  increases  from  the  origin  by  0.5°C  in  3 icm.  The 
response  of  this  layer  is  shown  in  Fig.  3«  The  com.bined  action  of 
increasing  surface  temperature  and  radiative  cooling  generates 
a fog  layer  which  grov/s  as  the  air  is  convected  over  the  viarm.er 
viater.  Near  the  origin  of  the  change  in  surface  tem.perature 
(x  = 100  m)  a thin  fog  layer  of  about  18  m depth  exists  with 
maximum  liquid  water  content  just  above  the  surface.  At  1000  m 
the  fog  layer  has  grown  to  about  40  m and  at  5000  m the  layer 
is  100  m deep.  At  this  depth  the  liquid  water  content  has 
developed  two  relative  maxima:  one  near  the  surface  and  one 
near  the  fog  top.  This  is  because  the  100  m deep  fog  bank 
is  no  longer  radlatively  thin  (Table  I)  and  fog-top  radiative 
cooling  is  beginning  to  dominate  over  surface  radiative  cooling. 
The  early  stages  of  a temperature  inversion  formation  at  the  top 
of  the  fog  bank  are  also  in  evidence. 

After  the  change  from  stable  to  unstable  conditions 
and  the  onset  of  radiative  cooling,  the  local 

shifts  to  a value  of  0.0008,  Such  a flux  ratio  places  a turbulent 
dominated  layer  in  the  v/arm  surface  fog-free  regime  of  Figure  1. 
Fog  nevertheless  exists  in  such  a layer  because  of  the  radiative 
cooling  present.  A convenient  measure  of  the  effect  of  radiation 
in  a fog  layer  of  thickness  5 is  an  effective  radiative  flux 


rad 


defined  as 


ELEVATION  ABOVE  SEA  LEVE 


Figure  3a.  Temperature  and  liquid  water  content  profiles  in  an 
advective-radiative  fog  at  x = 100  m . 
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LIQUID 


Figure  3b.  Temperature  and  liquid  water  content  profiles  in  an 
advective-radiative  fog  at  x = 1000  m . 


ELEVATION  ABOVE  SEA  LEVEL 


LIQUID 


Figure  3c.  Temperature  and  liquid  water  content  profiles  in  an 
advective-radiative  fog  at  x = 3000  m . 
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Direct  comDarison  of  0^  with  6 orovldesan  assessment 
of  the  strength  of  radiation  transport  compared  to  turbulent 
transport  in  the  fog  layer.  For  this  vrarrr.  surface  advection 
fog,  e is  approximately  0.5  . The  role  of  radiation 

is  therefore  pivotal  in  its  evolution. 

The  structure  of  radiative  cooling  in  the  fog  is  shown 
in  Fig.  4.  In  the  first  1000  m,  maximum,  radiative  cooling 
is  concentrated  at  the  surface  but  then  begins  to  shift  up- 
ward into  the  center  and  top  of  the  cloud  bank  as  it  becomes 
radiatively  thick  to  terrestial  radiation.  The  turbulence 
level  depicted  in  Fig.  5 shows  a steady  increase  in  x as 
the  unstable  mi.xed  layer  grows  with  the  maximum  t urbulence 
levels  concentrated  in  the  center  of  the  cloud  bank. 


r 
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6.  Boundary  layer  with  stratus  cloud 

A boundary  layer  over  a large  expanse  of  water  capped  by 
a subsidence  Inversion  will  turbulently  ingest  moisture  from 
the  surface.  Long  wavelength  cooling  of  the  v/ater  vapor 
laden  air  may  then  lead  to  condensation  somewhere  within  the 
layer.  Stratus  formation  brought  about  by  such  a sequence 
represents  one  of  the  simplest  and  most  basic  cloud  episodes 
possible  in  the  atm.ospherlc  boundary  layer.  In  an  expansive 
basin  or  region  influenced  only  by  advection  in  the  form  of  a 
steady  subsidence,  a periodic  steady  state  may  form  in  which 
a solar  driven  diurnal  cycle  exists  about  a mean  boundary  layer 
state  ’With  a steady  mean  inversion  height  controlled  by  the 
radiative  conditions  (season  and  latitude),  surface  temperature, 
and  upper  conditions  (wind,  water  content  and  subsidence)  set 
by  the  larger  circulation.  We  proceed  next  to  illustrate  such 
an  episode  and  to  relate  that  Illustration  to  occurrences  of 
expansive  stratus  layers  (v.'here  the  one-dimensional  unsteady 
version  of  turbulent -radiative  theory  would  be  applicable) 
such  as  the  summ^er  stratus  layer  along  and  off  the  coast  of 
California. 

The  conditions  for  this  illustration  are  selected  to 
generally  conform  to  the  conditions  for  sum^m.er  stratus  along 
the  California  coast  as  described  by  'lelberger  (19^^)  and 
more  recently  observed  by  Mack  et  al.  (197-).  Surface  bou-  ■ 
conditions  are  established  for  a fixed  tem.perature  se' 
of  17°C,  surface  water  m.ixing  ratio  at  saturati  , 
roughness  determ.ined  by  ^roude  scalinr  [Eo. 
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conditions  are  for  latitude  40°N  at  summer  solstice.  Geostrophic 
conditions  set  an  upper  level  wind  of  10  m/sec  and  a subsidence 
whose  characteristic  value  is  .5xl0"^sec"^  (Neiberger,  1961; 

Lilly,  1968).  Initial  conditions  (which  v/ill  be  lost  after 
several  days  of  simulation)  were  selected  to  correspond  to  a 

I 

mildly  stable  temperature  profile  of  lapse  30^/3z  = .003°C/m  j 

and  clear  sky  with  an  initial  relative  humidity  of  0.9.  The  | 

calculation  is  begun  shortly  before  sunset  the  first  day  and  | 

runs  for  six  days.  Behavior  during  day-light  hours  must  be 
considered  somewhat  approximate  in  the  present  illustration 
because  droplet  scattering  of  direct  solar  radiation  would 
reduce  solar  penetration  into  the  cloud  interior.  ■ 

The  quasi-perlodic  evolution  of  cloudiness  in  the  boundary 
layer  is  shown  in  Pig.  6.  The  stratus  grows  both  downward  and  I 

upward  during  nocturnal  periods  and  thickens  until  sunrise  when  j 

the  solar  heating  begins  to  disperse  the  stratus.  Radiative 
cooling  is  dominant  at  the  cloud  top  as  shown  in  Fig.  7,  and 
produces  an  unstable  lapse  rate  within  the  cloud  and  a corres- 
pondingly enhanced  turbulence  production  which  is  also  maximum 
in  nocturnal  periods  and  in  the  upper  portion  of  the  cloud 

I 

(Fig.  8).  The  radiation  cooling  at  the  cloud  top  strengthens  ; 

i 

the  capping  inversion  v:hich  oscillates  up  and  down  over  the  ! 

period  of  the  diurnal  cycle  (Fig.  9). 

We  note  that  the  predicted  cloud  base  in  these  diurnal 
cycles  is  highest  in  late  afternoon  and  lov;est  in  early  morning 
(Fig.  9)  while  cloud  top  is  highest  in  early  morning  and  lowest 
in  late  afternoon.  Correspondingly  the  inversion  height  is 
highest  in  early  morning  and  lov/est  in  later  afternoon — a 
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Figure  9.  Diurnal  variation  of  cloud  base  and  top  and  heat  and  moisture  fluxes  at  the  surface 
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result  which  is  opposite  to  that  for  inversion  height  cycles 
driven  by  solar  heating  in  the  cloud-free  boundary  layer 
(Lewellen,  et  al.,  197^).  This  predicted  result  is  in  accord 
with  the  stratus  observations  of  Nelberger  (19^^)  as  well  as 
those  of  Mack  et  al.  (197^)  in  which  they  detected  a diurnal 
variation  of  inversion  height  in  the  presence  of  California 
coastal  stratus  which  regularly  showed  maxima  in  the  early 
morning.  It  is,  of  course,  the  radiative-turbulent  drive  of 
cloud  top  cooling  which  allows  the  radiating  layer  to  propagate 
condensation  upward  during  the  night  while  turbulence  cools  the 
cloud  interior  below. 

Some  summary  aspects  of  this  illustration  in  comparison 
with  the  observations  of  Mack  et  al.  (197^)  which  guided  the 
choice  of  boundary  conditions  are  as  follows.  The  typical 
liquid  water  content  of  the  theoretically  generated  stratus 
was  0.1-0. 2 g/kg  compared  to  0.1-0. 3 g/kg  observed.  The 
average  height  of  stratus  base  was  50-200  m (theoretical) 
compared  to  0-200  m (observed).  The  average  height  of  stratus 
tops  at  maximum  inversion  were  600-800  m both  predicted  and 
observed.  The  theoretically  predicted  cooling  rate  100  m 
below  cloud  top  between  1 a.m.  and  3 a.m.  local  time  for  the 
last  two  days  of  the  simulation  was  0.023®C/hr  compared  to 
values  of  0.025®C/hr  between  1 a.m.  and  5 a.m.  observed  by  Mack, 
et  al.  The  amplitude  of  the  vertical  diurnal  motion  of  the  cloud 
base  was  100-150  m predicted  compared  to  observations  of 
50-150  m . 
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7.  Surface  fog  resulting  from  the  lowering  of  stratus 

We  have  seen  that  radiative-turbulent  coupling  of  sub- 
sidence capped  stratus  layers  propagates  the  base  downward 
during  nocturnal  cooling.  If  the  base  is  propagated  com- 
pletely to  the  ground,  a fog  is  formed  at  the  surface.  This 
class  of  fog  event  is  common  in  coastal  California  during 
periods  of  strong  subsidence  and  has  been  described  by 
Lelpper  (1948)  and  Mack,  et  al.  (1974).  Several  atmospheric 
fluid  mechanisms  have  been  postulated  to  explain  this  class 
of  fog  event  including  vertical  convection  driven  by  the 
sea  breeze-land  breeze  cycle  (Neiberger,  1944),  and  radiative 
cooling  of  stratus  during  the  night  (Mack,  et  al.,  1974). 

The  single  necessary  condition  all  observers  note  for  this 
class  of  fog  is  the  existence  of  the  late  afternoon  or  early 
evening  inversion  height  at  about  400  m or  lower.  The 
typical  Clear  summer  afternoon  temperature  profile  is  as 
indicated  in  Fig.  10,  as  measured  by  Mack  et  al.,(1974). 

A strong  upper  Inversion  exists  above  the  inversion  height 
z^  with  lapse  rate  typically  about  .04®C/m.  A transi- 

tion region  of  thickness  Az  = 40  m exists  above  z^  with 
a lapse  rate  of  .01®C/m  while  below  the  inversion  the  lapse 
rate  is  -.015®C/m.  We  begin  our  prediction  with  such  an 
initial  profile  over  a fixed  sea  surface  state  utilizing 
the  model  of  Parts  1 and  2 in  the  vertical  one-dimensional 
(z)  unsteady  version.  Initial  relative  humidity  is  .97  in 
the  region  up  to  the  inversion  falling  off  to  a value 


k ^ J 
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BtSr AVAILABLE  COPV 

corresponding  to  a mixing  ratio  of  q = .001  over  the  region 
of  the  transition  zone.  The  calculation  is  done  for  solar 
inclination  at  solstice  and  begun  at  7 p.m.  in  the  evening. 

Several  calculations  were  carried  out  with  the  height 
as  a parameter.  For  z^  = 300  m , the  atmosphere  begins  to 
cool  after  sunset  through  long  wavelength  radiation  brought 
about  by  the  water  vapor  in  the  air.  At  about  1 a.m.  the 
following  morning,  stratus  forms  extending  from  230  m to 
300  m . After  condensation  begins  the  stratus  grows  both 
downward  and  upward.  The  cloud  propagates  dovmward  to  within 

i 

about  l40  m of  the  surface  at  which  time  (5  a.m.)  the  cloud 
top  propagating  upward  exists  at  450  m . From  this  time  on, 
the  cloud  base  continues  to  lift  and  both  base  and  top 
propagate  upward. 

For  z^  = 200  m the  same  dynamic  occurs.  In  this  case,  | 

however,  the  downward  propagating  stratus  reaches  the  surface 
by  5 a.m.  and  remains  until  7 a.m.  when  the  cloud  base  lifts 
off  the  surface  (Fig.  11).  It  can  be  seen  from  Fig.  11 
that  there  is  an  Indication  that  surface  fog  forms  just 
before  the  stratus  base  reaches  the  ground  and  propagates 
upward  to  reach  the  dovmward  propagating  base,  the  two 
coalescing  together  at  about  10  m elevation.  During  this 
period  the  stratus  top  grows  from  230  m to  350  m . Through  the 
remainder  of  the  day  the  cloud  top  m.oves  upward  and  the  cloud 
base  lifts  until  by  mid-afternoon  the  top  exists  at  400  m while 
the  base  exists  at  225  m . 

As  can  be  seen  in  Fig.  12,  the  afternoon  (3  p.m.) 
cloud  has  thinned  (approximately  200  m deep)  compared  to  the 
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STRATUS  FORMATION  AND  LOWERING 


LIQUID  WATER  ISOPLETHS 


Figure  11.  Evolution  of  fog  resulting  from  the  lowering  of  stratus. 

Fog  shows  some  tendency  to  form  at  surface  just  before 
downward  propagating  stratus  reaches  the  surface. 


STRATUS  LOWERING  CASE 

LIQUID  CONTENT  PROFILE 


LIQUID,  G/KG 


Figure  12.  Variation  of  the  liquid  water  profile  in  a stratus  lowering 
fog.  Stratus  forms  at  about  1 am,  is  lowered  to  ground 
at  5 am,  and  lifted  at  3 pm  in  the  following  afternoon. 
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early  morning  (5  a.m.)  cloud  (approximately  350  m thick)  and 
the  average  liquid  water  content  is  reduced  about  50  percent. 

This  is  due  to  evaporation  brought  about  by  direct  solar 
’^adiation  during  the  daytime. 

The  measurements  of  .Mack  et  al.,  (1975)  :f  liquid  water 
content  in  stratus  lowered  to  the  surface  and  for  conditions 
typical  of  the  present  calculation  indicate  liquid  water  content 
ranging  from  0.1  to  0.2  g/kg  and  increasing  from  the  sur- 
face upward  into  the  cloud.  This  trend  is  consistent  with  the 
calculated  results  although  the  calculated  water  content  is 
less  than  measured  most  likely  due  to  the  presence  of  drizzle 
in  the  measured  case. 

In  Figs.  13  and  Ik  the  corresponding  distribution  of  radiant 
heating/cooling  is  shown.  Here  we  see  a manifestation  of 
the  large  difference  between  direct  solar  and  terrestial  ■ 
absorption  lengths.  Because  of  the  strong  terrestial  absorp- 
tion and  emission  of  the  water  droplets  the  cloud  top  is  in 
strong  radiative  cooling  v;hich  dominates  over  solar  heating 
even  in  daylight  hours.  Direct  solar  radiation  is  absorbed 
much  deeper  within  the  cloud  interior  and  this  deposited 
energy  transferred  as  heat  by  turbulence  to  the  surface. 
Evaporation  by  solar  heating  occurs  rather  uniformly  through- 
out the  cloud  because  of  the  long  absorption  length  for 
solar  radiation.  Strong  cloud  top  cooling  brought  about 


STRATUS  FORMATION  AND  LOWERING 


ISOPLETHS  OF  RADIATIVE  HEATING/COOLING 


Figure  13.  Radiative  cooling/heating  distribution  in  a stratus 

lowering  fog.  Maximum  cooling  is  concentrated  at  cloud 
top.  After  sunrise,  net  heating  occurs,  but  deep 
within  the  cloud. 
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STRATUS  LOWERING  CASE 
RADIATIVE  HEATING/COOLING  PROFILES 


QRAD=--^ 

/>Cp 


Figure  14,  Profiles  of  radiative  cooling/heating  at  1:30  am  (stratus 
just  forming) , 5 am  (stratus  on  the  surface) , and  3 pm 
(stratus  base  lifted) . 
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BEST  AVAILABLE  COPY 

by  the  short  terrestlal  absorption  length  provides  the 
Instability  and  turbulence  production  machinery  for  the  cloud 
to  turbulently  transport  the  heat  deposited  by  the  sun  In 
Its  Interior.  The  details  of  this  turbulent  transfer  can 
be  seen  In  Figure  15  where  the  turbulent  heat  flux  Is  shown 
directed  from  cloud  base  to  cloud  top.  It  should  also  be 
noted  that  heat  flows  turbulently  downward  above  the  cloud 
top  to  the  strong  radlatlvely  cooling  surface.  This  heat 
flux  Is  supported  by  the  turbulence  v/hlch  has  diffused  a 
few  tens  of  meters  above  the  cloud  top  Into  the  strongly 
stable  region  above  the  cloud. 

It  Is  noteworthy  that  the  corresponding  turbulent 
moisture  flux  shown  In  Figure  l6  Is  always  positive  representing 
a continuous  flow  of  water  from  the  sea  surface  upward  into 
the  cloud.  At  no  time  (after  stratus  formation)  during  the 
course  of  lowering  and  raising  was  the  water  flux  negative 
indicating  a diffusion  of  water  downward  from  the  cloud 
to  the  region  belov;  the  base.  Hence,  the  relative  downward 
flow  of  water  from  a radlatlvely  cooling  cloud  could  only 
take  place  by  precipitation  which  is  not  included  in  the 
present  calculation. 

Precipitation  transport  of  water  down  to  the  subcloud 
region  will  enhance  the  liquid  water  level  in  that  region.  The 
measured  liquid  water  content  shown  in  Fig.  12  is  greater  than 
predicted  most  likely  because  of  the  drizzle  present  in  the  mea- 
sured case.  The  present  theory  does  provide  some  insight  as  to 
the  appropriateness  of  the  entrainment  hypotheses  upon  which 
simple  cloud  top  mixing  models  are  based  (Lilly,  I966;  Schubert, 


0 
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1976).  In  Figure  15  the  distribution  of  the  turbulent  heat 
flux  in  diurnally  varying  stratus  Indicates  that  the  virtual 
heat  flux  is  positive  throughout  the  mixed  layer  and  has  a 

\ 

% 

positive  definite  average  over  the  layer  and  through  the 
cloud  top  where  a small  region  of  negative  flux  exists.  Thus, 
Ball's  hypothesis  that  the  average  should  vanish  (Bail,  196I) 
and  the  corresponding  maximum  entrainment  condition  v/hlch 
follows  from  it  is  not  descriptive  of  a cloudy  layer.  The 
present  turbulent- radiative  theory  indicates  that  although 
the  virtual  heat  flux  undergoes  significant  diurnal  variation 
(nearly  vanishing  when  the  sun  is  at  its  maximum  radiance),  it  is 
always  directed  upward  from  sea  to  cloud  in  the  sub-cloud 
region.  Further,  the  long  wavelength  radiative  heat  flux 
is  about  10  to  20  percent  of  the  turbulent  virtual  flux 
in  the  night  and  much  greater  than  the  turbulent  flux  during 
sunlight  hours  in  the  sub-cloud  region.  Hence,  on  average 
radiative  transport  is  as  important  as  turbulent  transport 
in  the  sub-cloud  region. 

The  overall  turbulence  levels  are  shown  in  Figure  17. 

In  the  early  morning  when  stratus  covers  only  a small  portion 
of  the  initial  mixed  layer  the  turbulence  levels  have  a typical 
dry  mixed  layer  structure.  After  the  cloud  base  rises,  the 
confinement  of  radiative  cooling  to  the  cloud  top  and  the 
presence  of  heating  in  the  cloud  interior  skevv's  the  distribu- 
tion to  concentrate  and  locate  the  maximum  turbulence  levels 
near  the  cloud  top. 

The  evolution  of  the  temperature  profiles  shov/n  in 
Figure  18  indicate  a steadily  intensifying  cloud-top  inver- 
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STRATUS  LOWERING  CASE 
TURBULENCE  PROFILE 


q2,  M^/SEC^ 


Figure  17.  Turbulence  levels  in  a stratus  lowering  fog. 
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STRATUS  LOWERING  CASE 
TEMPERATURE  PROFILE 


Figure  18.  Temperature  profile  evolution  in  a stratus  lowering  fog. 
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sion  compared  to  the  weaker  inversion  characteristic  of  the 
clear  sky  initia]  state.  An  inversion  of  approximately 
0.15°C/m  in  the  early  morning  before  sunrise  is  con- 
sistent with  the  measurements  of  Mack  et  al.  (197^). 

It  is  of  interest  to  follow  the  evolution  of  Q »/0 

* V* 

for  a stratus  lowering  fog  as  shown  in  Figure  19.  Since 
the  atmosphere  cools  over  a fixed  temperature  surface,  this 
fog  falls  into  the  warm  surface  class.  It  can  be  seen  that 
the  surface  is  Initially  far  form  critical.  During  the  night, 
(j^/6v^  falls  rapidly  towards  the  critical  value  but  never 
quite  reaches  it;  yet  surface  fog  does  exist.  Turbulent 
existence  criteria  do  not  apply  fully  here  because  of  the 
radiation  cooling  present.  If  vre  define  an  effective  0^,^  as 
by 

rad 


eff . 


'“v.> 


eff, 


= 0 +0^ 
V*  * 


,rad 


and  calculate  0''“'^  over  a region  5 where  the  heat  flux 

* 

is  approximately  constant,  we  obtain  a 15  percent  effect  of 

radiation  in  the  surface  layer  which  tends  to  elevate 

(q^/Q  ) . This  adjusted  value  is  denoted  (<?  *70  ) 

crlt.  ■ eff. 

It  can  be  seen  that  such  a radiative  adjustment  brings  the 

actual  Q into  the  surface  fog  regime  during  most 

of  the  period  for  which  the  cloud  base  was  at  the  surface. 

This  calculation  Indicates  that  radiation  cooling,  and 


the  corresponding  turbulence  production  it  effects  at  a 
stratus  top, cools  the  atmosphere  and  propagates  the  cloud  top 
upward  and  cloud  base  downward.  The  effective  region  over 
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STRATUS  LOWERING  CASE 


Odoy  I day 

Time,  days  8 hours 


Figure  19.  The  evolution  of  q*./6v*  ^ stratus  lowering  fog. 

Correction  of  (q.^/9y.^) " for  radiation  leads  to  (<?*/9v*)eff 
locating  the  lowered  stratus  in  the  warm-surface  fog 
regime  of  Fig.  1. 
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which  the  cloud  base  may  be  propagated  away  from  the  top  by 

I 

I 

turbulence  production  at  the  top  during  the  nocturnal  period 
is  determined  by  the  long  v/avelength  radiation  absorption 

i 

length  and  by  the  turbulent  scale  length  over  which  turbulence  ! 

may  be  diffused  downward  from  the  top  to  the  base.  For  j 

i 

these  typical  stratus  conditions  used  as  the  basis  of  the 
calculation,  we  find  the  region  to  be  350  m (Figure  13). 

The  empirical  observation  that  the  downv/ard  propagating 
stratus  base  can  reach  the  ground  during  nocturnal  cooling  if 
the  cloud  top  is  at  400  ra  or  lower  is  quite  consistent  with 
this  result. 
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8.  Concluding  Remarks 

The  interaction  between  turbulent  and  radiative  trans- 
port in  the  development  of  low  level  cloud  and  fog  has  been 
explored  using  a second-order  closure  model  of  turbulence. 

We  find  that  existence  criteria  for  turbulent  dominated  sur- 
face fog  may  be  established  on  the  basis  of  a critical  tur- 
bulent humitity  to  heat  flux  ratio.  Thermal  radiation  is  shown 
to  play  a significant  role  in  such  atmospheric  fog  and  stratus 
cloud  events.  Illustrations  of  the  dynamics  of  a turbulent, 
radiatively  coupled  lov;  level  atmosphere  for  warm  water  fog 
and  subsidence  capped  stratus  show  consistent  quantitative 
features  found  in  observations  of  such  fog  and  cloud.  Inversion 
heights  v/lth  cloud  in  the  boundary  layer  are  predicted  to  have 
maxima  in  early  morning  and  minima  in  late  afternoon  in 
contrast  to  the  diurnal  behavior  of  inversion  heights  in  a 
dry  boundary  layer.  Stratus  is  shovm  to  propagate  both  down- 
v/ard  and  upward  during  nocturnal  cooling  periods  under  the  in- 
fluence of  the  radiant/turbulent  driving  mechanism  located  in 
the  cloud  top.  If  the  late  afternoon  inversion  height  is 
low  enough,  the  stratus  is  predicted  to  reach  the  surface. 

The  important  effect  of  the  large  difference  between  solar 
and  terrestial  radiative  absorption  coefficients  for  liquid 
drops  on  the  turbulent  behavior  of  cloud  and  fog  is  made  mani- 
fest in  these  illustrations.  That  disparity  in  radiative 
Interaction  is  one  of  those  small  but  vital  differences 
characteristic  of  the  planet.  It  is  somev/hat  akin  to  the 
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density  maximum  of  water  at  ^®C  which  allows  ice  to  float  and 
thus  has  enormous  implications  for  bodies  of  water  on  the  earth' 
surface.  If  the  solar  and  terrestial  radiation  absorption 
lengths  bore  the  opposite  relationship  to  one  another,  cloud 
tops  would  be  heated  by  the  sun  while  the  cloud  interiors 
cooled  to  space  leading  to  stabilization  and  the  cessation  of 
turbulence  production  by  the  cloud  top  cooling  mechanism. 
Evaporation  of  the  cloud  interior  down  to  the  base  would  be 
delayed  until  radiative  "burning  away"  of  the  stratus  progressed 
downward  devoid  of  the  mixing  Influence  of  turbulent  transport. 
Cloud  structure  and  cloud  dynamics  would  be  significantly 
altered  and,  we  imagine,  so  too  would  our  atmosphere. 


s 
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APPENDIX  B 


COMPARISON  OF  THE  A.R.A.P.  SECOND-ORDER  CLOSURE  MODEL 
WITH  MEASUREMENTS  IN  A l-IARINE  BOUNDARY  LAYER' 


We  now  compare  the  predictions  of  the  A.R.A.P.  second-order 
closure  model  with  measurements  of  the  fluid  dynamic  behavior  of 
a marine  boundary  layer.  A carefully  studied  data  case  to  be  used 
as  the  basis  of  the  comparison  is  that  of  Pennell  and  LeMone  for 
the  trade-wind  boundary  layer  (Pennell  and  LeMone,  1974).  These 
measurements  were  taken  for  the  mean  profiles  of  velocity,  tempera- 
ture and  water  mixing  ratio  and  the  corresponding  turbulent  stress, 
heat  flux  and  moisture  flux  fields.  A full  comparison  of  the 
A.R.A.P.  predictions  of  the  dynamics  and  evolution  of  the  boundary 
layer  would  require  measurements  of  the  fluid  state  at  an  initial 
time  and  at  some  subsequent  time.  Unfortunately  such  sequential 
measurements  are  not  presently  available.  We  are  thus  only  able  to 
provide  a verification  of  the  consistency  of  the  A.R.A.P.  model 
in  the  sense  that  the  turbulence  co-existing  with  the  measured 
mean  profiles  is  well-represented  by  the  turbulence  correlations 
predicted  by  the  A.R.A.P.  model  and  based  upon  the  measured  mean 
profiles.  A more  complete  verification  of  the  A.R.A.P.  model 
must  await  the  availability  of  turbulence  and  mean  profile  data 
taken  over  several  different  time  levels. 

In  Figs.  1 to  3 we  indicate  the  mean  profiles  of  velocity, 
temperature,  and  total  water  mixing  ratio  taken  from  Pennell  and 
LeMone  (1974)  and  Sommeria-Klein  and  LeMone  (1977)  for  the  NCAR, 
Puerto  Rico  experiment.  The  data  points  were  used  as  initial  con- 
ditions for  the  A.R.A.P.  predictions.  The  sea  surface  temperature 
was  held  fixed  at  25.9°C.  The  reference  time  of  1200  hours  denotes 
that  the  solar  radiation  input  corresponds  to  starting  the  calcula- 
tion at  noon.  The  mean  profiles  evolve  relatively  slowly  until 
significant  cloudiness  sets  in  after  approximately  four  hours. 

The  consistency  of  the  model  predictions  can  be  checked  by  com- 
paring the  turbulent  flux  profiles  during  this  slow  evolution  of 
the  mean  profiles  with  those  measured  concurrently  with  the  mean 
profiles.  Figure  4 shows  the  variation  in  the  shear  stress  over 
a five  hour  simulation  period.  In  the  mixed  layer  the  predictions 
agree  quite  well  with  the  measurements.  Note,  in  particular,  that 
the  model  has  no  adjustable  input  at  the  surface.  The  effective 
surface  roughness  depends  internally  upon  the  predicted  uw  . 

Thus  the  agreement  of  the  surface  value  with  the  data  is  an  im- 
portant achievement  for  the  model.  There  is  some  discrepancy 
between  the  observations  and  predictions  at  the  top  of  the  mixed 
layer  where  an  undershoot  is  predicted  but  not  observed.  The 
' undershoot  intensifies  when  clouds  develop. 

; In  Fig.  5 we  show  the  same  shear  stress  profiles  but  also 

1 include  a predicted  result  for  this  case  based  on  the  NCAR  subgrid 

tThese  comparisons  were  made  by  Owen  Cote  while  on  leave  from 
j Staff  Meteorology,  Air  Force  Systems  Command,  Hamscom  AIB,  Bedford, 

Mass. 
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turbulence  model  (Sommeria-Klein  and  LeMone , 1977).  It  can  be 
seen  that  their  model  also  predicts  an  undershoot  which  incriases 
with  cloud  activity. 

In  Fig.  6 the  vertical  velocity  variance  is  shown.  There 
is  general  agreement  between  both  model  predictions  and  the  data 
in  the  mixed  layer  below  600  m altitude.  The  character  of  the 
A.R.A.P.  prediction  and  the  NCAR  prediction  above  the  mixed  layer 
is  similar  in  that  both  show  a strong  increase  in  a„  as  clouds 
develop . 

The  virtual  potential  temperature  flux  is  shown  in  Fig.  7. 
Again,  there  is  quantitative  agreement  between  both  predictions 
and  the  data  in  the  mixed  layer  and  qualitative  agreement  between 
the  two  predictions  above  the  mixed  layer. 

As  discussed  in  Section  3.1  the  biggest  discrepancy  between 
the  ITCAR  model  results  and  the  A.R.A.P.  results  is  the  prediction 
of  liquid  water  content  after  the  clouds  develop.  There  is  from 
1 to  2 orders  of  magnitude  difference  in  the  predicted  value  of 
average  liquid  water  content.  This  difference  remains  to  be 
resolved. 

Similar  calculations  for  comparison  with  data  from  GATE  day 
September  15  (LeMone,  1977)  have  also  been  made.  The  results 
of  this  comparison  were  quite  similar  to  the  ones  shown  here  for 
the  Puerto  Rico  case  but  with  somewhat  more  scatter  in  the  data. 

We  show  the  fluxes  of  momentum,  temperature  and  water  vapor  in 
Figs.  8 and  9. 

We  reiterate  that  these  comparisons  can  only  be  used  to  show 
a general  consistency  of  the  model.  The  comparisons  cannot  be 
used  to  make  a careful  assessment  of  the  accuracy  of  the  model. 
This  will  necessitate  a more  detailed  specification  of  the  experi- 
mentally-observed boundary  conditions. 
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Figure  1.  Evolved  mean  profiles  of  velocity  according  to  the  A.R.A.P 
model  and  the  model  of  Sommeria-Klein.  Initialization  of 
the  models  is  for  the  mean  velocity  measurements  shown  as 
Puerto  Rico  initial  input. 
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Figure  2.  Evolved  virtual  potential  temperature  profiles  predicted 
by  the  A.R.A.P.  model  for  the  conditions  of  Fig.  1. 
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Figure  3.  Typical  mean  profile  (total  water  mixing  ratio)  taken  from 
the  measurements  of  Pennell  and  LeMone  and  used  as  an  ini- 
tial condition  for  the  A.R.A.P.  model  (solid  line).  Sub- 
sequent evolution  of  the  mean  total  water  mixing  ratio  is 
shown  at  2 hrs  and  5 hrs  after  initialization  as  predicted 
by  the  A.R.A.P.  model.  After  5 hrs,  clouds  have  appeared 
in  the  modeled  boundary  layer. 
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Figure  5.  Comparison  of  measured  longitudinal  momentum  flux  with 

the  A.R.A.P.  model  predictions  and  the  subgrid  turbulence 
model  predictions  of  Sommeria-Klein  and  LeMone. 
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Figure  6.  Vertical  velocity  variance  predicted  by  the  A.R.A.P.  model, 
predicted  by  the  model  of  Sommeria-Klein,  and  the  NCAR 
1972  Puerto  Rico  measurements. 
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Figure  7.  Virtual  potential  temperature-velocity  correlations  as 
predicted  by  the  A.R.A.P.  model,  the  model  of  Sommeria- 
Klein,  and  compared  with  the  NCAR  1972  Puerto  Rico  mea- 
surements . 
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APPENDIX  C 


AIULYSIS  OF  A FRONT  MOVING  IN  THE 
PLANETARY  BOUNDARY  LAYER 

by 

H.  Segur 


The  occurrence  of  atmospheric  fronts  in  the  mid-latitudes 
is  a common  meteorological  event.  These  fronts,  which  carry 
with  them  relatively  sharp  temperature  discontinuities,  exist 
both  on  a synoptic  scale  (e.g.,  polar  fronts)  and  on  much 
smaller  scales  (e.g.,  thunderstorm  gust  fronts  and  sea  breeze 
fronts).  During  the  period  of  the  present  contract,  we  have 
begun  to  study  the  detailed  local  structure  of  an  atmospheric 
front,  using  a combination  of  analytical  and  computational 
methods.  In  this  section,  we  present  some  of  our  preliminary 
results . 

On  a cyclonic  scale,  the  development  of  warm  and  cold 
fronts  was  studied  numerically  by  Kasahara,  Isaacson  and  Stoker 
(1965) , using  a two-layer  model  of  an  incompressible  atmosphere 
on  a rotating  earth,  with  a B-plane  approximation.  Their  model 
had  two  (horizontal)  space  dimensions,  and  they  assiamed  hydro- 
static pressure  in  the  vertical  direction.  Further,  they  ignored 
the  dynamics  of  the  upper  layer.  The  result  was  a credible 
description  of  the  horizontal  motion  of  these  fronts,  on  a 
cyclonic  scale.  They  did  not  discuss  the  local  structure  of 
the  surface  front,  which  requires  a separate  analysis.  Some 
of  the  motivation  for  the  present  study  is  to  provide  a descrip- 
tion of  this  local  structure,  to  complement  their  description 
of  the  large-scale  dynamics. 

A numerical  study  of  the  fronts  associated  with  a thunder- 
storm gust  front  (Teske  and  Lewellen,  1977)  is  described  in 
Section  3.2  of  this  report.  A major  difference  between  these 
two  kinds  of  fronts  is  that  Coriolis  forces,  which  have  a 
dominant  effect  on  the  cyclonic  scale,  are  relatively  unimportant 
for  this  application.  In  either  case,  there  is  a need  for  a 
local  description  near  the  surface  front. 

A.  A Simplified  Analytical  Model  (without  Coriolis  forces) 

We  consider  the  two-dimensional  motion  of  two  layers  of 
incompressible  fluid,  moving  under  gravity  on  a horizontal 
surface,  as  shown  in  Fig.  1.  Here  h(y,t)  is  taken  as  the 
height  of  the  front  and  H(y,t)  may  be  taken  as  the  height  of 
the  tropopause.  To  make  this  first  model  as  simple  as  possible, 
we  neglect  Coriolis  forces.  This  creates  no  problems  for  the 
thunderstorm  gust  front,  but  does  limit  the  range  validity  of 
the  analysis  when  applied  to  a cold  front.  We  return  to  this 


point  below,  and  discuss  the  circumstances  in  which  these 
.esults  may  be  applied  to  the  motion  of  a cold  front.  The 
relevant  equations  of  motion  are  (in  Fluid  I) : 


3v  ^ 
3y  3z 


0 


lx  + 

3t 


lx  + i iE  = 

3z  p 3y  3z 


(1) 

(2) 


(in  Fluid  II) : 


lx' 

3t 


+ v' 


lx' 

3y 


lx'  4-  lx' 

3y  3z 


0 


+ w’ 


lx'  + 1 + !li 

3z  p (1-6)  3y  3z  3y 


(3) 

(4) 


where  t;  refers  to  the  shear  stress,  including  both  laminar 
and  turbulent  terms.  A fundamental  assumption  of  this  analysis 
is  that  the  pressure  is  purely  hydrostatic.  Thus,  we  assume 
the  pressure  is  constant  above  z = H(y,t)  ; in  fluid  II, 

p’  = gp(l-6)  (H(y,t)  - z)  + p^  : (5) 

and  in  Fluid  I, 

p = gp(l-6)  (H(y,t)  - h(y,t))+  gp(h(y,t)  - z)  + p^  .(6) 

The  validity  of  this  assumption  is  discussed  at  relevant  points 
below. 

The  vertical  velocity  vanishes  at  the  ground  (z  = 0) . At 
the  interface  (z  = h)  and  free  surface  (z  = H)  there  are 
kinematic  conditions, 

1^  + v'  1^  = w'  , on  z = H(y,t)  , (7) 


on  z = h(y, t)  . (8) 


We  anticipate  that  the  shear  stress  will  be  negligible  out- 
side of  this  boundary  layer  along  the  ground  and  near  the  inter- 
face, and  neglect  (t2,t0.  If  we  integrate  Eqs.  (1)  and  (2)  _ 
in  the  z-direction,  making  use  of  the  other  equations,  the  final 
result  may  be  written  as 


where  t.  is  the  stress  at  the  interface;  and  t is  that  at 
the  ground.  The  only  approximations  made  to  this^point  are  the 
hydrostatic  approximation  and  the  neglect  of  (9T/3y) . 

We  consider  first  the  case  in  which  we  may  neglect  shear 
stress  altogether.  Then  the  only  force  is  the  hydrostatic  pres- 
sure, which  introduces  no  z-variation  in  v (from  (2),  (4)). 
Thus,  3v/3z  vanishes  if  it  did  initially.  Equations  (9)  then 
reduce  to 

_ + _ (vh)  = 0 


i 

i 

i 


3v  > „ 3v  , 3H  , 3h  _ « 


+ 4 (v'(H-h))  - 0 


Iz'  4.  Iz' 

3t  3y 


For  atmospheric  applications,  6 is  relatively  small  (typically 
less  than  about  0.07). 


w 
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The  steady-state  version  of  this  model  was  used  by  Yih  and 
Guha  (1955)  in  a theoretical-experimental  study  of  internal 
hydraulic  jumps . They  found  good  agreement  when  6 was 
moderately  large  (v0.20).  For  6 = 0.05,  the  jumps  predicted 
were  about  30  percent  larger  than  those  observed;  no  more  than 
10  percent  of  this  discrepancy  could  be  attributed  to  scatter 
in  their  experimental  data.  It  appears  that  this  discrepancy 
is  a measure  of  the  validity  of  the  hydrostatic  pressure  ap- 
proximation. 

The  system  of  Eqs.  (10)  is  hyperbolic,  having  four  real 
characteristics.  Two  of  these  correspond  to  motion  in  which 
the  whole  system  responds  like  a homogeneous  fluid.  In  an 
atmosphere  without  shear  (v  = v'),  for  small  6 , these  speeds 
are  given  by 


^ = V ± /gE  (11) 

Further,  in  these  two  modes,  disturbances  satisfy 

A(v'-v)  =0  , A(h  - ^ H)  = 0 (12) 

so  that  the  density  discontinuity  at  z = h has  no  effect  to 
leading  order  in  6 . If  we  take  the  height  of  the  tropopause 
as  9 km,  propagation  speeds  for  this  mode  are  on  the  order  of 
200  m/sec,  much  faster  than  observed  speeds  of  atmospheric 
fronts . 

The  other  two  modes  of  the  system  (10)  correspond  to  waves 
traveling  along  z = h , and  these  are  of  more  interest  in  the 
present  study.  To  leading  order  in  6 (with  v = v'),  the 
characteristics  are  defined  by 


(13) 

and  disturbances  in  this  mode  satisfy 

AH  = 0 (14) 

Thus,  to  leading  order,  the  tropopause  does  not  feel  the  motion 
of  the  front,  and  the  upper  fluid  plays  a relatively  passive 
role.  In  a stationary  atmosphere  (v  = v'  =0)  with  a relative 
density  (or  temperature)  discontinuity  of  0.05  at  a height  of 
500  m,  wave  propagation  speeds  in  this  mode  are  on  the  order 
of  15  m/sec,  which  is  in  rough  agreement  with  observed  propaga- 
tion speeds  of  fronts.  Thus,  our  primary  interest  is  on  these 
internal  modes . 

The  model  of  Kasahara,  et  al . , 1965,  is  consistent  with 
(10)  except  that  they  also  included  (and  were  primarily  in- 
terested in)  Coriolis  effects.  We  will  see  below  that  (10)  is 
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a valid  approximation  to  their  model  for  phenomena  that  occur 
over  a few  minutes  (i.e.,  for  times  very  much  shorter  than  a 
rotation  period  of  the  earth) . They  also  assumed  that  the 
dynamics  of  the  upper  fluid  could  be  neglected.  Some  justifi- 
cation for  this  assumption  can  be  inferred  from  the  fact  that 
the  tropopause  responds  only  weakly  to  motion  of  the  front, 
according  to  (14) . This  additional  assumption  reduces  (10)  to 


II  + ^ (vh)  = 0 


(15) 


which  can  be  solved  by  characteristics,  and  which  corresponds 
exactly  to  the  eqviations  of  nonlinear  shallow  water  waves,  with 
a reduced  gravity  (g6) . 


Moreover,  in  this  simplified  model,  the  motion  of  a cold 
front  corresponds  to  the  motion  of  a wave  of  water  advancing  on 
a dry  bed  after  a dam-break  (cf..  Fig.  2).  This  problem  has 
a simple  solution,  given  by  Ritter  (1892);  for  -/g6ho  < y/t 
< 2/g6ho  , 


V = 2/3(/g6ho  + y/t) 
/g6h  = 1/3 (2/g6ho  - y/t)  , 
w = -2/3  z/t  , 


(16) 


Dressier  (1952)  noted  that  the  velocity  of  the  tip  of  the  for- 
ward wave,  2/g6ho  , exceeded  the  observed  velocities  by  as 
much  as  40  percent  (in  water-air  experiments,  where  5=1). 
He  argued  that  even  though  shear  stress  is  negligible  outside 
a "thin"  boundary  layer,  that  layer  engulfs  the  fluid  region 
near  the  tip,  where  h vanishes.  Thus,  the  shear  stress  at 
the  ground,  which  was  dropped  after  (9),  must  be  retained  to 
model  the  motion  of  the  front  accurately.  Dressier  (1952) 
and  Whitham  (1955)  used 


T = K v‘ 
e 


(17) 


where  K is  an  empirical  constant,  so  that  (15)  is  replaced  by 


3h  , J_ 
9t  3y 


(vh) 


= 0 


3v  , 
3t 


3v 

3y 


3y 


Kv^ 

"TT 


= 0 


(18) 


For  K <<  1 , this  additional  term  has  little  effect  except 
near  the  tip  (i.e.,  the  surface  front)  where  it  is  dominant. 


(a)  Initial  conditions,  t = 0 


— -yg8ho[— ayi 
(b)  Inviscid  solution,  t>0 


Figure  2.  Idealized  solution  of  the  motion  of  a cold  front 
after  sudden  release,  without  dissipation  or 
Coriolis  forces. 
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Dressier  and  Whitham  each  solved  (18) , by  somewhat  different 
methods;  their  final  results  are  nearly  identical.  The  de- 
velopment below  follows  closely  that  of  Whitham  (1955) , 

First,  we  find  the  local  shape  of  the  surface  front,  as 
predicted  by  (18).  Let  y = a(t)  denote  the  position  of  the 
tip  obtained  by  solving  (18) ; (recall  that  if  K * 0 , 
a(t)  = 2/g6ho ' t) . Near  the  tip, 

y ^ a(t)  , V a(t)  , h -»  0 

If  the  particle  acceleration,  3v/3t  + v 3v/3y  , remains 
finite,  then  the  dominant  balance  near  the  tip  must  be 

g6  II  + 0 . (19) 

Thus,  near  y = a(t)  , 

h(y,t)  ||  (4)^(a(t)  - y)^/2  (20) 

The  characteristic  length  in  this  region  is  the  radius  of 
curvature  at  the  tip : 


This  may  be  taken  as  a measure  of  the  height  of  the  front  near 
the  tip. 

Dressier  (1952)  had  previously  observed  from  available 
experimental  evidence  (for  the  dam-break  problem)  that  the  slope 
of  the  front  becomes  vertical  at  the  tip.  Equation  (20)  con- 
firms this  result  and  parameterizes  the  local  shape  in  terms 
of  the  density  (or  temperature)  difference,  6 , the  surface 
resistance,  K , and  the  front  speed,  a(t)  . This  prediction 
of  the  shape  of  the  surface  front  is  one  of  the  major  results 
of  this  model. 

It  might  be  appropriate  to  remind  the  reader  at  this  point 
that  (18) , as  well  as  whatever  results  may  be  derived  from  (18) , 
are  based  on  a hydrostatic  pressure  approximation,  and  that 
the  results  of  Yih  and  Guha  (1955)  suggest  that  this  approxima- 
tion loses  its  validity  as  6 ->■  0 . It  follows  that  the  current 
model  should  be  more  reliable  when  the  temperature  jump  is 
stronger  (provided  that  the  Coriolis  effects  can  still  be 
neglected) . 

Even  for  strong  fronts,  however,  vertical  accelerations 
must  be  important  wherever  the  slope  of  the  front  becomes 
vertical,  as  predicted  at  the  tip  by  (20).  Thus,  the  hydro- 
static pressure  approximation  (obtained  l5y  neglecting  vertical 
accelerations)  must  break  down  near  the  tip  of  the  front,  for 
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any  6 . However,  this  objection  applies  with  equal  validity 
to  the  dam-break  problem  (6  = 1) . The  success  of  the  theory 
in  that  application  suggests  that  the  effect  of  this  breakdown 
may  not  be  too  severe. 

Next,  we  need  to  determine  the  tip  velocity,  a(t)  . 

Whitham  (1955)  argued  that  (16)  applies  except  in  a region  near 
the  tip,  where  h is  too  small  to  neglect  the  resistance  term 
in  (18),  (cf..  Fig.  3).  Within  this  dissipative  region, 

(C(t)  < y < a(t)),  the  velocity  is  assumed  to  be  constant  (in  y) , 
The  mass  flux  into  this  region  is 

M = ^ = phj^(V-0  ; (22) 


where  hi(t)  is  the  height  of  the  front  at  y = 5(t) 
we  will  match  to  (16) . From  (16)  , 

1=2  V(t)  - /g6ho 

C = 2 V - /pHo  + I Vt  • 

h,  (t)  = h (l ^ 

^ ° y 2/g6ho  j 


where 


(23) 


After  substituting  (23)  into  (22)  and  integrating,  one  obtains 
the  total  mass  in  the  dissipative  region: 

3 

M(t)  = h/g^  ( I - t . (24) 

° V 2v^g<5ho ' 

The  momentum  flux  across  y = C(t)  being  pVh]^(V-5)  , 
the  rate  of  change  of  momentum  in  this  region  is  given  by 

^ = pVhj^(V-C)  + forces  . (25) 

a)  Since  we  are  neglecting  the  dynamics  of  the  warm  air 

mass,  the  net  horizontal  pressure  gradient  acting  on  this 
region  is  99 

r ph^  p(l-6)h^  1 

s - g 2 — J • 

(The  caveats  noted  above  about  the  hydrostatic  pressure  approxi- 
mation  apply  here  as  well.) 

b)  The  drag  force  along  the  ground  is 

a 

pK ^ v^dy  = pKV^(a-5) 
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c)  We  have 


dP 


^ (MV)  = VM  + IW 


d)  At  the  tip,  V = a , and  this  must  apply  throughout 
this  region.  Substituting  all  of  these  into  (25)  and  rearranging 
yields 


1 

1 


• -7^  t = 

2*'g6hoy  /g6ho 

K/g6ho/  g N ^ r a 

\/g6ho/  [_/g6ho 


1 - 


2/g6ho 


- t 


(I 


a 


/g6ho 


•■>] 


(26) 


We  note  that  the  inviscid  solution, 
a(t)  = 2/g6Tr  t 

satisfies  (26)  if  K = 0 . For  K # 0 , it  determines  how  the 
motion  of  the  surface  front  is  modified  by  shear  stress  at  the 
ground.  The  equation  may  be  put  in  dimensionless  form  by  de- 
fining a dimensionless  time. 


T 


K 


(27) 


and  a dimensionless  lag  of  the  tip  behind  its  inviscid  position, 


(2/g6h^  t - a)K 


Then  (26)  becomes 

4(a)^  xd  + = 16 (2-a) ^ (3aT-2a) 

and  the  position  of  the  tip  is  given  by 

h 

a = 2/g6h^  ^ “ X ‘ 


(28) 

(29) 

(30) 


We  are  primarily  interested  in  the  solution  of  (29)  for 
relatively  short  times  (on  the  x-scale) . Whitham  (1955)  showed 
that 

a(T)  = 2.58916  x^^^  [1  + 0.81917  x^^^  + ©(x^^^) j . (31) 

In  terms  of  our  original  variables,  the  position  and  velocity 
of  the  surface  front  are  given  by 
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a(t) 


a(t)  * 2/g6h 


2/g^  t - 2.5816  (g6h^)^/^  * 

jjL  + 0.81917  + . . . J 

o[l  - 1-726U 

- 1.76747  + ...  J 


(32) 


(33) 


while  the  position  of  the  rear  of  the  dissipative  region  is 
given  by 


C(t) 


/g6H^  ^ 

- 5.3024  + . . . J 


(34) 


These  trajectories  are  plotted  in  Fig.  4 for  a set  of  typical 
parameter  values . 


This  completes  the  development  of  this  (simplest)  model. 
Next,  we  consider  some  of  its  consequences.  First,  we  note  that, 
according  to  this  simple  theory,  the  motion  of  a cold  front 
moving  under  a stationary  warm  air  mass  in  the  absence  of 
Coriolis  forces  is  essentially  unsteady.  As  sho^m  in  Fig.  4, 
the  surface  front  continually  decelerates  from  its  initial 
velocity.  Whitham  (1955)  noted  that 

a V t^^^  , a t~^^^  as  t -►  » (35) 


so  that  the  front  neither  comes  to  rest  nor  attains  a steady 
velocity.  Apparently  the  dynamic  effects  of  the  warm  air 
mass  (including  Coriolis  effects)  are  required  to  obtain  a 
stationary  front  (or  a warm  front) , as  we  show  below. 


Second,  let  us  consider  the  effects  of  Coriolis  forces, 
which  we  have  neglected  to  this  point.  If  we  still  retain 
9/9x  = 0 , which  is  legitimate  near  the  center  of  a cold  or 
warm  front,  (18)  should  be  replaced  by 


9h  , 

9t  9y 


(vh) 


9v 

9t 


+ V 


9h  . Kv^  _ 


(36) 
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where  u is  the  velocity  component  parallel  to  the  front. 
(There  is  an  additional  equation  for  3u/3t  , which  we  omit 
here.  Moreover,  if 


u >>  V , 

the  resistance  term  also  should  be  modified  Co  act  on  the  total 
kinetic  energy.)  Under  what  circumstances  can  the  Coriolis 
term  in  (36)  be  neglected? 

The  time-scale  of  the  inviscid  rotation-free  solution 
(K  = 0 , f = 0 in  (36))  is 


t^^  'v-  /fTTI^  (37) 

The  time-scale  of  the  resistance  term  in  (36)  is 


^2 

h/Kv  , 

which  depends 

on  the  local  values  of  v and  h . 

Near  the 

tip. 

g'S 

V = a 'V  /g6h^ 

» 

from  (21).  Thus, 

(38) 

It  should  not 

be  surprising 

that  these  time-scales 

align,  since 

the  resistance  term  exactly  balances  the  pressure  gradient  near 
the  tip.  The  relevant  time-scale  for  the  Coriolis  term  in  (36) 
is 

t^  '''  I fu/v]  (39) 

To  neglect  this  term,  it  is  necessary  that 

ti  -v  t2  <<  t^  . (40) 

For  a fairly  weak  front. 


6 .01 


and  under  any  reasonable  circumstances, 

h < 9 X 10^  m 
o 


Thus , 


t2  < 300  seconds  (41) 

in  the  worst  case,  and  100  seconds  is  probably  more  typical. 
Taking 
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1/f  10  seconds  , 

we  see  that  (40)  requires  that 

lu/v|  <<  30  . (42) 

The  numerical  results  of  Kasahara,  et  al.  certainly  satisfy 
(42) , and  it  appears  that  (42)  may  be  satisfied  by  many  cold 
fronts . 

These  estimates  should  not  be  interpreted  to  mean  that  the 
current  model  can  be  applied  directly  to  cold  fronts  generated 
by  Coriolis  forces.  Rather,  they  suggest  that  the  dynamics 
predicted  by  this  model  near  the  tip  typically  occur  much  faster 
than  the  time-scale  on  which  Coriolis  forces  act.  Thus,  the 
behavior  near  the  tip  should  be  relatively  insensitive  to 
Coriolis  effects,  provided  hydrostatic  forces  dominate  Coriolis 
forces  away  from  the  tip. 


B.  Stationary  Fronts 

We  noted  above  that  our  simplest  model  (without  Coriolis 
forces)  produced  only  cold  fronts,  in  which  the  cold  air  con- 
tinually advanced  beneath  the  warm  air.  Moreover,  the  surface 
front  was  always  blunt-nosed,  as  shown  in  Fig.  5.  If  Coriolis 
forces  are  included,  warm  fronts,  cold  fronts  and  stationary 
fronts  are  all  possible,  as  demonstrated  by  Kasahara,  et  al . 
(1965).  In  this  section,  we  analyze  the  structure  at  the  sur- 
face  of  a stationary’’  front,  for  comparison  with  that  in  a cold 
front.  We  will  see  that  the  blunt-nosed  shape  that  occurred 
above  cannot  occur  in  a stationary  front,  but  that  the  surface 
front  may  either  be  needle-nosed  or  have  a finite  slope. 


If  we  include  Coriolis  forces  (with  a 0-plane  approximation), 
Eqs.  (10)  become 


I?  ^ - It  ^ If  + = 0 


3x 


If  ^ If  + i - 0 


I?'  + u’  I?'  + V'  II  + g II  + fu'  = 0 , 


3t 


3x 


(43) 


|H'  + u'  |H'+  V'  |H'+  g |§ 


3C 


9x 


3y 


3x 


fv'  = 0 


^ (H-h)  + ^ (u'(H-h))  + ^ (v'(H-h)) 


(c)  stationary  front  with  Coriolis  forces.  Worm 
air  in  motion 


Figure  5.  Possible  shapes  of  surface  fronts. 


1 
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where  f = 2 oj  sin  (p  is  the  (constant)  Coriolis  parameter, 


exact 

solution : 

V 

= 

C 

II 

o 

u’ 

= 

const . 

u 

= 

const. 

3H 

fu' 

3y 

g 

3h 

1 

1 — 1 

p 

1 

3y 

(44) 


If  X is  in  the  east-west  direction,  this  solution  represents 
a stationary  polar  front.  [It  is  known  to  be  unstable  to 
cyclonic  perturbations,  which  grow  on  a time-scale  of  hours. 
However,  our  primary  interest  is  with  the  structure  of  the 
surface  front;  we  have  seen  that  the  relevant  time- scale  there 
is  on  the  order  of  a few  minutes.  Therefore,  we  may  ignore 
phenomena  that  take  much  longer  to  occur.] 

The  solution,  (44),  is  subject  to  the  same  criticism  that 
applied  to  (16):  near  the  surface  front,  where  h ^ 0 , the 
stress  at  the  ground  must  affect  the  dynamics  significantly. 
Thus,  the  first  two  equations  should  be  replaced  by 


3v 


3v 


^ + u ^ + V 


3t 


I?  + « + V + g(l-S)  |a  + g6  |£  - fv  . - -2#^ 


3x 


3u 


3v 

3y 

3v 


+ g(i-a)  |2  + g6  «l+fu  = 


3y 


3H 


3y 


3h 


3t 


3x 


3y 


3x 


3x 


(45) 


Our  objective  here  is  to  find  a suitable  modification  of  (44) 
which  approximately  satisfies  (45)  , even  as  h -*■  0 . 

An  immediate  generalization  of  (44)  is  obtained  by  letting 
(u,u’)  depend  on  y . Thus,  there  is  a family  of  solutions 
of  (43) , defined  by 

V = v'  = 0 , 

u’  = u'  (y) 
u = u(y) 

= - IhI 

9y  ~ g 

= - X [u 

3y  g6  ^ 


(46) 


(l-6)u'] 
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The  previous  solution,  (44),  is  a member  of  this  family.  We 
want  to  choose  u(y)  and  u' (y)  so  that  (46)  remains  a valid 
approximate  solution  of  (45) , even  as  h -*■  0 . 

Because  v'  = v h 0 in  (46),  we  must  have 

Ty,  = 0 . (47) 

For  consistency  with  (17) , we  define 


T = K u u 

xz,g  g'  ' 


Tx2,i  = -Kilu+u' I (u-u')  . (48) 


Ordinarily,  K > 


We  consider  first  the  case  in  which  the  warm  fluid  is  at 


rest  : 


u'  (y)  = 0 


= 0 
3y 


In  order  for  the  stress  terms  on  the  right  side  of  (45b)  to 
remain  bounded  as  h vanishes , we  require  that  as  h -►  0 , 


u V at 


p > 1/2  . 


Substituting  this  into  (46)  yields 


dh  „ fa  , p 

3y  • i6  ^ 

We  obtain  a valid  solution  only  if 

1/2  < p < 1 


h V C(y-y^)-‘-"P  , H = const.. 


u V D(y-y  ) ^“P  , u'  = 0 


In  this  case,  the  slope  of  the  surface  front  is  balanced  by 
the  (vanishing)  velocity  difference  across  the  front. 

To  this  order  of  approximation,  the  exponent  p is 
restricted  only  to  a finite  range.  For  any  p in  this  range, 
however,  the  surface  front  has  an  elongated  "needle-nose," 
as  shown  in  Fig.  5.  This  shape  contrasts  sharply  with  the 
blunt-nosed  shape  obtained  above  when  Coriolis  forces  were 
neglected  (also  shown  in  Fig.  5). 

Next,  suppose  that,  in  the  vicinity  (in  y)  of  the  surface 
front, 

u'(y)  Cl'  f 0 . (52) 


J 
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Now  the  stress  terms  in  (45b)  are  bounded  only  if 


K. 


^1/2 


u u 


K.+K 
V X g> 


as  h 


(53) 


Here  the  local  solution  is 


_ f 
3y  g 


ill  = X 

3y  ~ g6 


(1-5)  - 


K. 

1 


a/2-1 


^Ki+K 


(54) 


This  configuration  is  shown  in  Fig.  5,  as  well.  Again,  the 
slope  of  the  surface  front  is  balanced  by  the  velocity  dif- 
ference across  the  jump.  In  this  case,  the  velocity  difference 
does  not  vanish  near  the  surface  front,  and  the  slope  of  the 
front  remains  finite.  The  slope  away  from  the  surface  front 
can  be  either  larger  or  smaller  than  that  given  by  (54) , 
depending  on  the  local  velocity  difference.  For  sufficiently 
small  5 , the  slope  of  the  surface  front  is  opposed  to  and 
much  larger  than  the  slope  of  the  tropopause  above  it.  This 
feature  seems  to  be  frequently  observed  in  atmospheric  fronts. 


We  may  summarize  our  results  about  the  possible  shapes  of 
surface  fronts  as  follows: 


(a)  Fronts  accompanying  events  whose  duration  is  sig- 
nificantly less  than  an  hour  or  so  can  be  analyzed  without 
considering  Coriolis  effects.  These  fronts  are  always  cold 
fronts.  The  shape  of  the  surface  front  is  blunt-nosed,  as 
determined  by  a balance  between  the  hydrostatic  pressure 
gradient  and  the  resistance  of  the  ground  to  the  motion  of  the 
front  itself. 


(b)  For  events  which  last  for  several  hours,  Coriolis 
forces  can  play  an  important  role.  Warm  fronts,  cold  fronts 
and  stationary  fronts  are  all  possible. 

(c)  The  slope  of  a stationary  front  at  the  surface  is 
balanced  by  the  local  velocity  difference  across  the  front, 
subject  to  the  constraint  that  there  be  no  unbalanced  shear 
stress  from  motion  parallel  to  the  front.  The  slope  of  the 
front  at  the  surface  is  proportional  to  the  vertically-averaged 
velocity  in  the  warm  air.  This  is  also  true  at  the  tropopause. 
If  the  velocity  parallel  to  the  front  in  the  warm  air  does  not 
vanish,  the  slope  of  the  surface  front  is  opposed  to  and  usually 
larger  than  that  of  the  tropopause  above  it.  If  it  does 
vanish,  so  that  the  warm  air  is  at  rest,  the  slope  of  the  sur- 
face front  vanishes,  yielding  a needle-nosed  shape . The  sur- 
face front  cannot  be  blunt-nosed,  because  there  is  no  shear 
stress  normal  to  the  front. 


C . Comparison  with  Observations 


In  light  of  the  preliminary  and  somewhat  restricted  nature 
of  these  results,  an  exhaustive  comparison  of  the  predictions 
of  these  models  with  the  observed  behavior  of  surface  fronts 
is  thought  to  be  premature  at  this  time.  Nevertheless,  it  may 
be  worthwhile  to  show  that  some  of  the  qualitative  features 
predicted  by  these  models  are  observed,  even  when  our  somewhat 
restrictive  assumptions  (either  no  Coriolis  forces  or  stationary 
fronts)  are  not  satisfied.  This  qualitative  agreement,  in 
circumstances  beyond  the  range  of  the  present  model,  suggests 
that  the  models  can  be  generalized  to  a wider  class  of  problems 
without  affecting  the  qualitative  features  discussed  above. 

Figure  6 shows  three  vertical  cross-sections  of  the  at- 
mosphere over  the  United  States,  taken  from  an  analysis  of 
Hurricane  Hazel  by  Palmen  (1958) . These  figures  each  were  taken 
at  a fixed  latitude  (either  40°N  or  45®N),  with  east  to  the 
right  of  the  figure.  The  vertical  scale  is  given  on  the  figure 
but  the  horizontal  scale  is  not;  consequently  it  is  difficult 
to  measure  frontal  slopes  for  any  quantitative  comparisons. 

As  noted  in  the  captions,  the  heavy  solid  lines  denote  fronts 
or  tropopauses,  and  these  are  our  main  concern.  In  Figs.  6a 
and  6b,  the  cold  air  mass  is  in  the  center  of  the  figure,  with 
warm  air  on  either  side;  in  Fig.  6c,  the  cold  air  is  to  the 
left.  Temperature  differences  across  the  front  were  on  the 
order  of  - 10°C  (6  .02  - .03).  Five  surface  fronts  are 

seen  in  these  figures;  each  is  moving  to  the  east  (as  inferred 
from  the  wind  indicators).  Thus,  none  of  the  fronts  are  ac- 
tually stationary;  moreover,  Coriolis  forces  are  important  in 
determining  the  overall  motion  shown.  Thus,  neither  of  the 
models  developed  above  can  be  expected  to  apply  precisely. 

Even  so,  some  of  the  qualitative  features  of  the  fronts  are 
predicted  by  the  models,  as  we  now  show. 

Because  the  fronts  are  all  moving  eastward  (i.e.,  to  the 
right),  the  fronts  to  the  right  in  Figs.  6a  and  6b  are  cold 
fronts,  while  those  to  the  left  are  warm  fronts.  The  front 
in  Fig.  6c  is  a cold  front.  It  is  easily  seen  that  the  cold 
front  in  Fi^.  6a  is  blunt-nosed,  and  that  neither  of  the  warm 
fronts  in  Figs.  6a  and  b are.  (In  the  figures,  all  of  the  cold 
fronts  are  drawn  blunt-nosed  very  close  to  the  ground.  However, 
the  data  is  availdble  only  along  the  vertical  lines  shown,  and 
it  is  not  clear  that  the  shapes  of  the  cold  fronts  at  the  ground, 
as  shown  in  Figs.  6b  and  6c,  actually  were  observed.)  In  our 
model,  this  shape  is  a consequence  of  the  shear  stress  at  the 
ground,  which  slows  the  motion  of  the  cold  front.  The  cold 
front  in  Fig.  6b  is  almost  stationary  so  that  this  shear 
stress  is  relatively  small,  and  its  blunt  nose  should  be  ex- 
pected to  be  small  as  well.  Similarly,  the  speed  of  the  cold 
front  in  Fig.  6a  is  larger,  and  its  blunt  nose  ought  to  be 
larger  as  well.  Thus,  there  is  some  evidence  that  the  observed 
shapes  of  the  cold  fronts  at  the  surface  (blunt-nosed)  are 


Fir  ly  Vcriiol  cioM-»ec«i*>n,  OcUtbef  is.  »J«o  OCT,  ipprojtinutcly  along  pualkl  40*  N Thin  daihcd  Hnci 
are  >su(h«rmi,  and  heavy  linct  nuik  frontal  boundatiei  or  tropopauiei.  The  ihin  »olid  Unci  mark  itoiachi  at 
iniervali  of  ao  kl  of  obierved  wind,  without  regard  to  direction. 
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consistent  with  the  model,  as  is  the  dependence  of  the  nose 
size  with  front  speed.  No  attempt  at  quantitative  comparisons 
is  made  here. 

The  model  predicts  that  the  slope  of  a stationary  sur- 
face front,  the  slope  of  the  tropopause  above  it,  and  the 
vertically-averaged  velocity  (parallel  to  the  front)  of  the 
warm  air  between  are  all  related.  In  particular,  the  slope 
of  the  surface  front  is  opposed  to  the  slope  of  the  tropopause, 
and  typically  larger  than  it.  If  the  slope  of  the  tropopause 
vanishes,  the  slope  of  the  (stationary)  surface  front  beneath 
it  also  vanishes,  giving  the  surface  front  an  elongated  nose. 

The  opposing  slopes  of  the  front  and  the  tropopause  are  evident 
for  the  fronts  in  Fig.  6,  even  though  the  fronts  are  not  station- 
ary; (the  point  is  questionable  for  the  cold  front  in  Fig.  6a, 
where  the  slope  of  the  tropopause  is  changing  rapidly) . The 
tropopause  becomes  nearly  horizontal  west  of  the  warm  front 
in  Fig.  6b;  the  warm  front  beneath  it  is  seen  to  have  a rather 
elongated  tip,  as  predicted  by  the  model.  Thus,  the  qualitative 
behavior  of  these  fronts  is  in  agreement  with  the  model  of 
stationary  fronts. 

This  comparison  of  model  predictions  with  observed  atmos- 
perhic  fronts  seems  to  be  favorable,  at  least  in  terms  of 
predicting  qualitative  behavior.  Certainly  it  does  not  "prove" 
the  model.  It  does  suggest  that  the  model  has  some  predictive 
power,  and  should  be  able  to  provide  some  guidance  to  the 
detailed  computation  of  the  structure  of  fronts  in  the  planetary 
boundary  layer. 
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